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ABSTRACT:  We  present  a  method  for  calculating  elastic  and  inelastic  scattering  prob- 
abilities  for  light  particles,  such  as  helium  and  molecular  hydrogen,  scattering  from  surfaces 
with  which  they  weakly  interact.  The  method  is  a  unitary  one-phonon  approximation  in 
which  the  scattering  probabilities  are  calculated  from  thermally  averaged  amplitudes  which 
are  generated  numerically.  The  thermal  averaging  procedure  is  more  general  than  this  ap¬ 
plication  and  could  be  applied  to  other  systems  with  weak  inelastic  scattering.  We  also 
discuss  an  approximation  for  the  gas-surface  interaction  potential  that  can  greatly  simplify 
calculations  where  it  is  applicable.  Finally  we  present  some  preliminary  results  using  this 
method  to  study  rotationally  mediated  selective  adsorption  resonances  in  HD  scattering 
from  copper. 
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I.  INTRODUCTION 


In  all  experiments  in  which  helium  or  molecular  hydrogen  are  scattered  from  single 
crystal  surfaces  both  inelastic  scattering  and  selective  adsorption  resonances  are  present 
and  can  have  an  important  effect  on  the  elastic  scattering.  As  the  resolution  of  these 
experimemts  has  improved  not  only  is  it  possible  to  extract,  from  the  elastic  scattering, 
kinematical  information  like  diffraction  peak  positions  and  selective  adsorption  energies  but 
it  has  become  possible  to  extract  dynamic  information  like  the  diffraction  peak  intensities 
and  the  resonance  lineshapes.  It  also  has  become  possible  to  investigate  inelastic  scattering 
processes  leading  to  an  understanding  of  the  full  gas-surface  interaction.  To  interpret  these 
experiments  and  to  compare  the  results  with  calculated  potential  energies  it  is  necessary 
to  be  able  to  accurately  describe  the  dynamics  of  the  scattering  process.  In  this  paper 
we  present  a  method  for  calculating  scattering  probabilities  that  can  treat  both  inelastic 
scattering  and  selective  adsorption  resonances  in  the  presence  of  the  other.  By  studying 
inelastic  molecular  hydrogen  scattering,  in  particular  the  effects  of  the  rotational  degrees  of 
freedom  on  the  scattering  process,  we  hope  to  both  describe  these  scattering  experiments 
and  discuss  qualitatively  how  selective  adsorption  resonances  and  inelastic  scattering  affect 
each  other  in  more  general  scattering  situations. 

This  method  is  a  unitary  one-phonon  approximation.  Unitarity,  i.e.,  that  the  sum 
of  all  the  calculated  scattering  probabilities  is  one,  is  necessary  for  studying  the  thermal 
attenuation  of  the  elastic  scattering  probability  and  for  studying  the  coupling  of  selective 
adsorption  resonaces  and  inelastic  scattering.  As  we  showed  in  a  previous  paper1  inelas¬ 
tic  scattering  is  greatly  enhanced  by  selective  adsorption  resonaces,  so  much  so  that  the 
distorted  wave  Born  approximation  breaks  down.  This  breakdown  is  caused  by  the  over- 
counting  of  scattering  events  that  is  inherent  in  the  Born  approximation,  and  is  corrected 
for  in  this  new  approximation  by  removing  the  overcounting.  In  spite  of  the  breakdown  of 


the  distorted  wave  Born  approximation  at  selective  adsorption  resonances  a  one-phonon 
approximation  is  still  the  appropriate  approach  for  studying  helium  and  molecular  hydro¬ 
gen  scattering  because  away  from  resonance  the  inelastic  scattering  probability  is  still  weak, 
so  that  once  a  resonant  particle  scatters  inelastically  its  subsequent  scattering  probability 
is  low. 

In  the  experiments2  we  wish  to  describe,  the  scattering  of  H2,3  HD,4  and  He,5  the 
low  masses  and  moments  of  inertia  of  these  particles  cause  those  particles  that  diffract 
or  undergo  rotational  transitions  to  leave  the  surface  in  well  separated  directions.  This 
angular  separation  allows  the  detection  of  these  effects  simply  by  changing  the  relative  angle 
of  the  detector  with  respect  to  the  substrate  and  source.  Further,  since  the  physisorption 
potentials  wells  are  shallow,  the  bound  states  for  these  particles  are  well  separated  in 
energy  and  can  be  observed  in  these  experiments  through  selective  adsorption  resonances. 
Selective  adsorption  resonances  are  observed  in  the  intensity  of  outgoing  scattering  peaks 
as  a  function  of  the  incidence  condtions.  They  are  due  to  virtual  diffractive  or  rotational 
transitions  into  bound  states  and  are  caused  by  either  the  corrugation  of  the  surface  or  its 
translational-  rotational  coupling. 

While  the  inelastic  scattering  for  these  systems  is  weak  enough  to  allow  the  obseravtion 
of  coherent  elastic  scattering  effects,  it  is  not  negligible  nor  unobservable.  The  kinematical 
constraints  of  the  scattering  process  permit  single  Rayleigh  phonon  creation  and  absorp¬ 
tion  events  to  be  seen  in  the  time  of  flight  spectra  of  helium  atoms.5  Rayleigh  phonons 
are  normal  modes  of  a  semi-infinite  surface  that  are  localized  to  the  surface;  for  a  given 
wavevector  parallel  to  the  surface  they  have  an  energy  that  is  lower  than  all  the  bulk 
phonon  energies.  The  dispersion  relation  of  the  Rayleigh  mode  with  respect  to  the  par¬ 
allel  momentum  allows  the  identification  of  various  time-of-flight  peaks  as  single  Rayleigh 
phonon  transitions.  In  the  case  of  molecular  hydrogen  scattering  the  decreased  energy  res¬ 
olution  and  the  increased  inelastic  scattering  probability  make  it  more  difficult  to  observe 
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single-phonon  peaks  experimentally.  The  stronger  van  der  Waals  attraction  felt  by  the 
hydrogen  molecule,  compared  to  the  helium  atom,  causes  the  molecule  to  scatter  from  a 
steeper  part  of  the  repulsive  potential6  which  increases  the  inelastic  scattering.  Increased 
multiphonon  scattering  masks  single-phonon  peaks  in  the  time  of  flight  distributions.  An¬ 
other  aspect  of  inelastic  scattering  of  current  experimental  interest  is  the  degree  to  which 
the  rotational  degrees  of  freedom  of  molecular  hydrogen  affect  the  rate  of  sticking  in  either 
physisorption7  or  chemisorption®  states.  One-phonon  calculations  can  describe  the  initial 
trapping  step,  when  the  sticking  procedes  through  such  a  step.  But  to  completely  describe 
the  sticking  process  it  is  necessary  to  describe  the  transition  of  a  trapped  particle  into  a 
stuck  particle,9  and  that  is  beyond  the  scope  of  the  current  method. 

The  simple  scattering  behavior  of  helium  and  molecular  hydrogen  allow  straight¬ 
forward  calculation  of  the  scattering  rates.  The  small  observed  number  of  elastic  and 
internally  inelastic  channels  makes  coupled  channels  calculations  practical  while  the  low 
inelastic  scattering  probabilities  mean  one-phonon-change  calculations  are  relevant  to  the 
description  of  the  scattering  process.  For  heavier  atoms  and  molecules  both  the  number 
of  elastic  channels,  which  are  not  well  separated  experimentally,  that  have  to  be  included 
and  the  increased  inelastic  scattering  make  both  the  observation  and  the  calculation  of 
elastic  and  inelastic  scattering  more  difficult. 

In  many  ways  the  present  calculation  is  related  to  a  long  line  of  previous  calculations 
but  it  has  some  new  features.  To  describe  the  translational-rotational  coupling  we  use  a 
coupled  channels  approach  in  which  the  wavefunction  is  expanded  in  the  spherical  hamonics 
so  that  the  rotational  behavior  is  calculated  using  a  discrete  set  of  states  instead  of  a 
continuum.  The  coupled  channels  description  of  scattering  from  a  static  substrate  has  been 
used  to  very  accurately  describe  both  helium19  and  molecular  hydrogen11  experiments. 
The  main  point  of  such  calculations  has  been  to  extract  the  particle-surface  potential 
from  the  r  ■rimenta!  scattering  intensities.  Since  direct  inversion  of  experimental  data 
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to  generate  the  potential  is  impossible,  it  is  necessary  to  search  through  possible  potentials 
to  find  the  one  for  which  the  calculated  scattering  probabilities  best  match  those  observed 
experimemtally.  Alternatively  the  bound  state  energies,  extracted  from  slective  adsorption 
resonances,  can  be  used  to  determine  the  potential.  Provided  both  that  the  coupling  that 
leads  to  the  resonances  is  weak  and  the  resonances  are  well  separated  from  each  other, 
the  bound  state  energies  can  be  extracted  from  the  energies  of  the  resonances  using  the 
kinematics  of  the  scattering  process.  If  these  conditions  are  not  satisfied  then  either  a  full 
coupled  channels  calculation  or  its  equivalent  is  required  to  match  the  experimental  and 
calculated  resonces. 

Inelastic  gas-surface  scattering  probabilities  have  been  calculated  in  many  ways  rang¬ 
ing  from  classical12  to  fully  quantum  mechanical.  Between  these  two  extremes,  the  lat¬ 
ter  of  which  is  discussed  in  more  detail  below,  calculational  approaches  have  included 
wavepacket  calculations13  (both  for  elastic  and  inelastic  scattering),  and  various  semiclas- 
sical  calculations14  in  which  the  motion  of  the  scattering  particle  is  treated  classically  and 
the  phonons  are  treated  quantum  mechanically. 

For  helium  and  molecular  hydrogen,  quantum  mechanical  calculations  are  the  most 
valid  approaches  due  to  both  the  quantum  mechanical  nature  of  the  elastic  scattering  pro¬ 
cess  and  the  discrete  nature  of  the  one-phonon-change  scattering  process.  The  simplest 
quantum  mechanical  approach  is  the  distorted  wave  Born  approximation,13  which  is  a 
Fermi’s  golden  rule  approach.  The  rate  for  transitions  caused  by  the  creation  or  absorp¬ 
tion  of  single  phonons,  is  calculated  between  scattering  state  eigenfunctions  of  the  fiat 
static  surface.  The  eigenfunctions  that  are  used  in  these  calculations  depend  only  trivially 
on  all  of  the  degrees  of  freedom  except  the  motion  normal  to  the  surface.  This  golden 
rule  approach  can  be  used  to  include  higher  order  phonon  processes  but  unless  further  ap¬ 
proximations  are  made  the  resulting  calculations  are  difficult  to  carry  out.16  As  computing 
power  has  increased  distorted  wave  Born  approximation  calculations  have  been  carried  out 


using  eigenfunctions  in  which  several  of  the  degrees  of  freedom  are  coupled,  for  instance 
eigenfunctions  of  corrugated  static  surfaces17  and  surfaces  with  translational-rotational 
coupling.1  The  approach  we  are  using  in  this  paper  is  related  to  these  extended  distorted 
wave  Born  approximation  calculations  in  that  we  treat  the  rotational  diffraction  on  an 
equal  footing  with  the  elastic  scattering. 

The  distorted  wave  Born  approximation  is  less  appropriate  for  calculating  the  elastic 
scattering  probabilities  when  they  are  changed  by  the  possibility  of  scattering  inelastically. 
The  simplest  calculation  of  the  elastic  scattering  probability  under  these  conditions  is 
to  multiply  the  scattering  probabilities  by  a  Debye- Waller-like  factor.18  Another  relatively 
simple  calculation  is  to  add  a  phenomenological  local  optical  potential19  to  the  rigid  surface 
potential  when  doing  a  coupled  channels  calculation.  The  optical  potential  simulates  the 
transfer  of  scattering  probability  from  the  elastic  scattering  to  the  inelastic  scattering  by 
absorbing  intensity  from  the  elastic  scattering  probability.  This  optical  potential  aproach 
can  be  improved  by  solving  for  the  non-local  energy-dependent  self-energy20  that  correctly 
describes  the  elastic  scattering  probability.  A  related  approach  is  to  include  inelastic 
scattering  in  a  scattering  matrix  calculation  to  study  its  effect  on  the  elastic  scattering 
lineshapes.21 

A  unitary  calculational  scheme,  in  which  all  of  the  scattering  intensity  is  accounted  for 
and  sums  to  unity,  involves  an  extended  coupled  channels  calculation,22  extended  in  the 
sense  that  possible  inelastic  transitions  are  included  in  the  coupled  channels  calculation. 
The  calculation  is  done  for  several  initial  occupations  of  the  lattice  modes  and  then  the 
scattering  probabilities  arc  thermally  averaged  with  respect  to  these  initial  occupations  to 
get  scattering  probabilities  to  compare  with  experiment.  Since  the  phonon  modes  form 
continua  and  any  number  of  phonons  can  be  created  or  destroyed,  it  is  necessary  to  truncate 
the  set  of  inelastic  processes  that  are  allowed. 

The  present  calculation  is  related  to  both  the  optical  potential  approach  and  the  ex- 


6 


tended  coupled  channels  approach  and  forms  a  bridge  between  the  two.  It  differs  from 
the  extended  coupled  channels  approach  in  that  the  thermal  averaging  is  done  before 
the  wavefunctions  are  calculated;  it  also  differs  in  that  the  solutions  are  calculated  itera¬ 
tively  instead  of  at  one  time.  The  final  results  of  this  calculation  for  the  elastic  scattering 
probability  are  identical  to  those  of  a  self-energy  calculation20  mentioned  above  but  are 
calculated  using  different  intermediate  quantities.  In  this  calculation  we  chose  to  truncate 
the  allowed  phonon-changes  to  include  only  one-phonon-change  processes,  but  include  all 
the  phonon  modes  in  the  lattice.  We  also  specify  how  to  treat  the  situation  in  which  par¬ 
ticles  can  trap  on  the  surface,  a  situation  that  can  cause  great  difficulty  if  sufficient  care 
is  not  taken. 

We  apply  this  method  to  HD  scattering  for  two  main  reasons,  the  most  important 
being  that  experiments  have  been  performed  on  this  system.  The  other  reason  is  that 
with  present  computer  power  this  method  is  most  applicable  to  HD  scattering.  Significant 
savings  in  computation  time  result  from  being  able  to  ignore  the  corrugation  of  the  surface, 
as  is  discussed  in  Sec.V.,  and  from  treating  broad  as  opposed  to  narrow  selective  adsorption 
resonances,  because  the  iterative  calculations  converge  faster.  To  study  selective  adsorption 
resonaces  and  inelastic  scattering  simultaneously  and  to  take  advantage  of  the  time  savings 
mentioned  above  requires  studying  HD  scattering. 

This  paper  is  organized  as  follows:  section  II  gives  the  general  derivation  of  the  ther¬ 
mal  averaging,  section  III  contains  a  unitary  one-phonon-change  approximation  using  the 
results  of  the  previous  section,  section  IV  discusses  the  use  of  stationary  state  scattering 
results  to  calculate  the  scattering  probabilities,  section  V  gives  the  approximations  that 
made  on  the  form  of  the  gas-surface  interaction  potential,  section  VI  presents  the  result  of 
some  preliminary  calculations  using  this  method,  and  section  VII  contains  a  summary  of 
the  main  results  of  this  paper. 

n 

t 


II.  THERMAL  AVERAGING 


In  this  section  we  show  how  thermally  averaged  scattering  probabilies  can  be  calculated 
in  terms  of  amplitudes  that  are  already  thermally  averaged.  This  procedure  does  not 
violate  our  concepts  of  statistical  mechanics  because  the  amplitudes  for  the  particle  are 
not  thermally  averaged  with  respect  to  the  particle  but  with  respect  to  the  phonons  from 
which  the  particle  is  scattering.  We  are  averaging  a  reduced  time  evolution  operator 
with  respect  to  the  phonon  coordinates;  the  time  evolution  operator  has  been  reduced  by 
operating  it  on  the  initial  state  of  the  particle.  We  call  this  reduced  time  evolution  operator 
an  amplitude-operator  throughout  this  text  to  remind  the  reader  that  it  is  both  an  operator 
on  the  phonon  coordinates  and  an  amplitude  for  the  scattering  particle.  The  thermally 
averaged  amplitudes  are  the  set  of  the  thermal  averages  of  products  of  this  amplitude- 
operator  with  all  possible  combinations  of  phonon  creation  and  destruction  operators.  This 
amplitude-operator  can  be  used  to  calculate  all  the  properties  of  the  particle  scattering 
from  any  surface  so  that  all  the  thermally  averaged  properties  can  be  calculated  from  the 
thermally  averaged  amplitudes.  These  thermally  averaged  amplitudes  obey  a  hierarchical 
set  of  equations  of  motion  which  can  be  solved  and  used  to  give  the  thermally  averaged 
scattering  without  explicitly  thermally  averaging. 

Our  development  of  the  amplitude-operator  is  based  on  studying  the  time  evolution 
of  the  state  of  the  system.  We  start  at  some  initial  time  (/  =  0)  with  the  scattering 
particle  in  some  initial  state  that  is  localized  sufficiently  far  from  the  surface  so  that  it 
is  not  interacting  with  it,  then  the  state  of  the  system  is  given  by  a  product  ket  of  that 
initial  state,  |x(0)),  and  the  initial  state  of  the  phonons, |{n,}),  which  is  specified  by  the 
occupation  of  all  the  normal  modes  of  the  surface 

l*(0))  =  lx(0)>|{n,}>.  (2.1) 

The  straight-forward  method  of  calculating  the  thermally  averaged  scattering  probabilities 


is  to  calculate  the  scattering  probabilities  for  the  initial  particle  state  scattering  from  each 
lattice  otate  in  an  ensemble  of  surfaces  and  then  averaging  the  scattering  probabilities 
weighted  by  the  thermal  probability  of  each  surface  in  the  ensemble.  Instead  of  this 
procedure,  we  show  that  we  can  calculate  the  thermally  averaged  scattering  directly  from 
thermally  averaged  amplitudes. 

The  time  evolution  of  the  state  of  the  whole  system,  from  which  the  scattering  prob¬ 
abilities  can  obtained,  is  given  by  operating  on  the  state  (2.1)  with  the  exponetial  of  the 
full  Hamiltonian  multiplied  by  the  time 

l*(0>  =  e'**‘!#(o).>  (2.2) 

Since  the  Hamiltonian  couples  the  motion  of  the  scattering  particle  with  that  of  the 
phonons  the  state  is  no  longer  be  a  product  state  once  the  particle  starts  to  interact 
with  the  surface  but  can  be  thought  of  as  a  sum  of  product  states.  Below  we  write  this 
sum  of  states  in  terms  of  how  the  occupations  of  the  phonon  modes  have  changed. 


A.  HAMILTONIAN. 

The  Hamiltonian  is  broken  into  three  terms:  the  Hamiltonian,  of  the  lattice  in  the 
absence  of  the  scattering  particle,  the  Hamiltonian,  of  the  particle  in  the  absence  of 

coupling  to  the  phonons,  and  the  potential, Ttn(,  that  couples  the  motion  of  the  scattering 
particle  with  the  motion  of  the  phonon  coordinates.  The  particle  Hamiltonian  includes 
both  the  kinetic  energy  of  the  particle  and  the  potential  that  couples  the  particle  with  the 
static  lattice;  the  interaction  term  is  the  potential  that  couples  the  particle  with  the  lattice 
minus  the  interaction  with  the  static  lattice 

H  =  Hlat  +  Hrart  4-  Vtnt.  (2.3) 
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The  weak  inelastic  scattering  and  small  changes  in  phonon  occupations  allows  us  to  describe 
the  lattice  by  a  harmonic  Hamiltonian  because  the  subsequent  anharmonic  effects  are  even 
smaller  in  a  real  system.  In 

Hlat  =  X^WAfllaA'  (2-‘4) 

A 

the  phonon  modes  are  labelled  by  a  composite  index  A  which  includes  the  momentum  of 
the  mode  parellel  to  the  surface,  the  polarization,  and  either  the  asymptotic  momentum 
normal  to  the  surface  far  from  the  surface  for  bulk  derived  modes,  or  the  decay  length  into 
the  surface  for  modes  that  are  localized  to  the  surface. 


In  this  section  we  make  several  approximations  only  to  simplify  the  presentation.  For 
example,  we  use  a  Hamiltonian  for  the  scattering  particle  that  just  consists  of  the  kinetic 
energy  due  to  its  center  of  mass  motion  and  a  potetial  energy  due  to  the  presence  of  the 
static  lattice  that  depends  on  the  position  of  the  center  of  mass 

2 

Hpart  =  ^  +  V(r).  (2.5) 

The  momentum  of  the  center  of  mass  is  p  and  the  postion  of  the  center  of  mass  is  r. 
Throughout  this  paper  lower  case  bold  face  letters  refer  to  three-dimensional  vectors  (e.g. 
r),  bold  face  upper  case  letters  refer  to  two-dimensional  vectors  that  are  perpendicular  to 
the  surface  normal  which  we  choose  to  be  the  z-direction  (e.g.  R),  and  italic  versions  of 
the  bold  face  letters  are  the  norm  of  the  vector  (e.g.  p). 


For  simplicity  we  take  only  that  term  in  the  interaction,  Vrtnf,  between  the  particle  and 
the  phonons  that  is  linear  in  the  phonon  coordinates 


= 


—  V 

^  A 


Vx{r)a[  +  h.c. 


(2.6) 


The  factor  of  jY-1/2  associated  with  the  sum  over  the  modes  of  the  lattice  (N  is  the 
number  of  atoms  in  the  surface)  comes  from  writing  the  displacement  of  each  lattice  atom 
in  terms  of  the  normal  coordinates  of  the  lattice  (see  section  5  Fq.  (5. 4-5. 6)  ). 
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The  procedure  for  thermal  averaging  that  we  present  does  not  depend  on  either  the 
neglect  of  the  internal  degrees  of  freedom  of  the  scattering  particle  or  the  exclusion  of 
more  than  the  linear  term  in  the  interaction  potential.  We  make  these  approximations 
to  improve  the  clarity  of  the  presentation.  In  particular,  later  we  include  the  rotational 
degrees  of  freedom  to  study  molecular  hydrogen  scattering. 


B.  TIME  EVOLUTION. 

Since  we  are  not  interested  in  the  behavior  of  the  lattice  except  to  the  extent  it  affects 
the  motion  of  the  particle  we  define  a  quantity  'I'  that  we  refer  to  as  an  amplitude-operator. 
It  is  an  amplitude  for  the  particle  and  an  operator  for  the  phonon  coordinates  with  the 
time  dependence  of  the  lattice  factored  out 

^(r,0  =  <r|  e~tHt  |x(0))  .  (2.7) 

In  a  related  gas-surface  scattering  calculation  Celli  and  Maradudin20  use  this  amplitude- 
operator  to  calculate  elastic  scattering  probabilities  in  the  presence  of  inelastic  scattering. 
The  amplitude-operator  and  the  time  evolution  operator  of  the  uncoupled  lattice  give  the 
time  dependence  of  the  molecule  scattering  from  any  initial  set  of  phonon  coordinates  when 
they  operate  on  that  initial  state.  In  particular  using  Eqs.  (2.1,  2.2,  2.7),  (rl'I'(f))  can  be 
written  as 

<r|*(0)  =  e-‘//'“‘t^(r,0  !{«,})•  (2-S) 

The  utility  of  this  amplitude-operator  is  that  its  time  dependence  is  independent  of  the 
initial  state  of  the  lattice;  this  independence  allows  the  calculation  of  the  the  scattering 
from  any  particular  surface  in  the  ensemble  using  just  this  one  operator.  The  equation  of 
motion  for  'k, 

,  jU(r,()  -  J/parf  (r)'i'(r,  t)  +  \\nt(rj)*(r  ,t).  (2.9) 
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is  obtained  by  taking  the  partial  derivative  with  respect  to  time  of  the  defining  equation  Eq. 
(2.7).  The  two  terms  containing  the  lattice  Hamiltonian  cancel.  The  particle  Hamiltonian 
commutes  with  the  lattice  Hamiltonian.  The  time-dependence  of  the  interaction  potential 
is  that  due  to  the  non-interacting  lattice 


Wr.O  = 


-t  Htattv 
e  vtnte 


Vx(r  )elulxla[  +  h.c. 


\  L 


(2.10) 


The  time  dependence  of  Vint  is  not  that  of  any  of  the  conventional  approaches  (Schrodinger, 
Heisenberg,  or  interaction  picture)  but  is  the  interaction  picture  for  the  lattice,  a  picture 
in  which  operators  have  the  time  dependence  of  only  part  of  zeroth  order  Hamiltonian, 
i.e.,  that  of  the  lattice  Hamiltonian.  We  choose  to  treat  the  time  dependence  this  way 
because  we  are  only  interested  in  the  lattice  time  dependence  as  far  as  it  affects  the 
particle’s  motion.  Since  the  time  dependence  is  known  exactly  it  can  be  removed  from  the 
calculation. 


Since  the  lattice  Hamiltonian  commutes  with  any  particle  operator  and  does  not  change 
any  particle  state  when  operating  on  it,  we  can  use  the  amplitude-operator  to  evaluate  the 
expectation  value  of  any  particle  operator.  For  example,  consider 

<*(O|0Partl#(O)  =  /  d*r  f  dV  (K}|^t(r,t)Opart(r,r,)4'(r,,t)|{nJ),  (2.11) 

where  we  have  inserted  two  complete  sets  of  particle  states  (in  the  position  representation) 
around  the  particle  operator.  Note  that  Opar{(r,r')  can  be  taken  outside  of  the  phonon 
matrix  element.  Since  Eq.  (2.11)  is  valid  for  any  operator  Opar(,  all  the  properties  of  the 
scattering  particle  are  contained  in  the  density  matrix,  which  is  'f/  times  its  Hermitian 
conjugate.  Accordingly,  thermally  averaging  the  expectation  of  any  particle  operator  just 
involves  integrating  over  the  thermally  averaged  density  matrix 


{W)\Cp*rtm)))th=  I  j  d'r  Opart[ r,  r’)  (  ^(r,  t)  ^(r',  t))^  ■  (2-12) 


Calculating  the  thermally  averaged  density  matrix  from  quantities  that  are  independent 
of  the  initial  state  of  the  lattice  is  possible  because  the  time  evolution  of  the  amplitude- 
operator  is  constructed  to  be  independent  of  the  initial  state  of  the  lattice. 


C.  MULTI-PHONON  EXPANSION. 


To  further  simplify  the  calculation  of  the  scattering  probabilities  we  break  up  the 
amplitude-operator  into  amplitude-operators  for  each  of  which  there  is  a  specific  change  in 
the  occupation  of  the  lattice  modes.  In  particular  there  is  a  term  in  which  the  occupation 
of  each  of  the  lattice  modes  has  had  no  net  change,  a  sum  of  terms  in  which  exactly 
one  mode  has  had  a  net  increase  or  decrease  of  one  phonon,  and  further  terms  for  every 
possible  net  change.  We  write  these  terms  by  explicitly  factoring  out  all  of  the  unpaired 
phonon  operators  from  each  term.  Each  phonon  creation  or  destruction  operator  in  the 
interaction  potential  which  lead  to  these  changes  in  phonon  occupations,  has  associated 
with  it  both  a  sum  over  the  modes  of  the  lattice  and  a  factor  of  JV-1/2,  where  N  is  the 
number  of  atoms  in  the  lattice.  We  write  this  factor  explicitly  in  front  of  the  sum  so  that 
all  of  the  newly  defined  amplitude-operators  are  not  proportional  to  any  power  of  the  size 
of  the  system  (see  below).  Thus,  the  decomposition  of  the  amplitude-operator  in  terms  of 
phonon-change  amplitude-operators  is  given  by 

=*Oph{r,t) 


£ 


(0  ^  lph(r  ~h)  d"  aA  (0  ^  lph(r’  ) 


(2.13) 


Each  of  these  new  phonon-change  amplitude-operators  can  be  written  as  a  sum  of  terms 
that  are  amplitudes  (i.e.  they  contain  no  phonon  operators)  times  pairs  of  phonon  opera¬ 
tors  such  that  the  index  of  both  the  creation  and  the  annilation  operator  in  each  pair  is 
the  same.  Associated  with  each  pair  of  phonon  operators  there  is  a  factor  of  /V-1  and  a 
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sum  over  the  lattice  modes.  There  is  only  one  sum  with  each  pair  of  operators  because  we 
have  taken  only  those  terms  in  which  the  operators  are  paired;  the  rest  of  the  terms  from 
the  same  applications  of  the  interaction  potential  are  in  higher  order  phonon-change,  terms 
because  they  have  more  unpaired  phonon  operators.  For  example,  the  zero-  phonon-change 
term  can  be  written 

*0  Ph(r,0  =  V'S Pfc(r.0  +  Jj52'l>Oph(r't)a\aX  +  (214) 

A  A,A' 

All  of  the  amplitudes  and  amplitude-operators  discussed  in  this  paper  are  independent  of 
the  size  of  the  system  in  the  limit  that  the  size  of  the  system  becomes  infinite  because  all 
sums  of  paired  operators  over  3 N  terms  are  cancelled  by  a  factor  of  N~l.  Some  of  the 
sums  over  lattice  modes  of  unpaired  operators  should  be  restricted  so  as  not  to  include  the 
terms  that  are  counted  elsewhere  because  of  pairing,  but  the  corrections  are  not  important 
or  even  relevant  because  all  the  restricted  terms  are  negligible  in  the  large  N  limit. 


Since  any  expectation  value  over  phonon  states  requires  that  all  the  phonon  operators 
must  be  paired  to  give  a  non-zero  result,  the  thermally  averaged  density  matrix  can  be 
written  as  a  sum  of  the  thermal  averages  of  the  square  of  each  phonon-change  term 

(*t(r,()*(r',t))tk  =  {*Jp)l(r,l)*op).(r'.<)}lA 


th 


+  {^U(r^.^-)aAaA^lPA(r'^.A,-)) 


th 


+ 


N 2 


A, A' 


(2.15) 


All  of  the  creation  and  destruction  operators  on  the  right  hand  side  of  this  expression  are 
explicitly  paired,  and  all  are  summed  over.  In  each  of  these  expectation  values  there  are 
many  terms  that  contribute;  contributing  terms  are  called  correlated  if  the  index  of  one 
of  the  paired  phonon  operators  is  identicle  to  that  of  another  pair.  In  the  large  N  limit 
all  correlated  terms  have  a  factor  of  N~^  that  is  not  associated  with  a  sum  and  hence 
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vanish  as  the  number  of  lattice  atoms  becomes  infinite  (see  Appendix  A).  This  vanishing 
of  correlated  contributions  allows  each  of  these  thermal  expectation  values  of  a  product  of 
operators  to  be  written  as  a  product  of  expectation  values 

(**(r,<Wr'.<))M  =  <)),„  (*o,»(r'.f))tk 


th 


(2.16) 


A,V 

Although  the  left  hand  side  of  Eq.  (2.15)  and  the  first  term  on  the  right  hand  side  of  the 
same  equation  appear  similar,  the  difference  is  that  in  the  quantity  on  the  right  hand  side 
all  the  phonon  operators  in  each  factor  are  explicitly  paired  while  in  the  quantity  on  the  left 
hand  side  each  factor  contains  unpaired  operators  that  when  paired  with  operators  from 
the  other  factor  give  rise  to  the  rest  of  the  terms  on  the  right  hand  side  of  the  equation. 


D.  THERMALLY  AVERAGED  AMPLITUDES. 

Thermally  averaged  n-phonon-change  amplitudes  can  be  defined  by  the  thermal  expec¬ 
tation  value  of  the  n-phonon-change  amplitude-operators,  in  which  we  thermally  average 
the  paired  phonon  operators  in  the  particle  amplitude-operator.  These  thermally  averaged 
amplitudes  can  be  equivalently  defined  by  taking  the  thermal  expectation  of  the  product 
of  some  phonon  creation  and  destruction  operators,  which  have  the  time  dependence  of 
the  uncoupled  lattice,  with  the  amplitude-operator.  This  equivalence  results  from  there 
being  only  one  term  in  Eq.  (2.13)  that  does  not  have  any  unpaired  phonon  operators  when 
multiplied  by  the  product  of  phonon  operators.  The  zero-  phonon-change  amplitude 

tf0p/,(r,0  =  {*0Ph{r,t))lh  =  (2-17) 
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is  the  thermally  averaged  amplitude  for  the  particle  given  that  the  lattice  is  in  its  initial 
state.  We  show  that  this  amplitude  describes  the  elastic  scattering  of  the  particle.  There 
are  two  types  of  one-phonon  change  amplitudes,  those  that  are  the  thermally  averaged 
amplitudes  for  the  particle  given  that  the  occupation  of  exactly  one  phonon  mode  has 
increased  by  one 


0lPh(r,*,A,+)  3  (*lph(r,t,\,+))th  =  («a(0*( r.*))|h. 


(«A  +  0 


and  those  given  that  the  occupation  of  exactly  one  mode  has  decreased  by  one 


(2.18) 


<!> WMA -)  2  (*lpfc(r,i,  A, («j,(0*(r.<))(t  •  (2-19) 


These  amplitudes  and  the  further  n-phonon  change  amplitudes  describe  the  inelastic  scat¬ 
tering  of  the  particle.  The  density  matrix  can  be  written  in  terms  of  these  thermally 
averaged  amplitudes  and  the  thermally  averaged  occupation  numbers  of  the  lattice  modes 
as  (see  Eq.  (2.16)) 

(®j(r.  0(r't  t))tfc  =^oP/t(r’  Ph(r'’  0 

d"  )nX^’lph.ir  >  A'  — )  I 


(2.20) 


+ 


—  V.... 

/v2  ^ 


A, A' 


The  advantage  of  writing  the  density  matrix  in  terms  of  the  thermally  averaged  amplitudes 
is  that  these  amplitudes  obey  a  hierarchical  set  of  equations  of  motion  and  do  not  need  to 
be  explicitly  thermally  averaged. 

The  equation  of  motion  for  the  zero-phonon-change  amplitude,  0Op/i-  'n  which  none  of 
the  lattice  mode  occupations  have  changed,  is  found  by  thermally  averaging  the  equation 
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of  motion  (2.9)  for  the  amplitude-operator 


.  d 

t  — 

dt 


0Op/»M)  =  HP^t{r)rp0ph{r,t)  +  ^  ^ 


(nA  +  l)V,A(r)0lpfc(r,t,A,+) 


+  *\V\  {r)^iph{r,t,Xt  -) 


(2.21) 


The  equations  of  motion  for  the  one-phonon-change  amplitudes,  in  which  only  one  phonon 
mode  occupation  has  changed,  are  found  by  thermally  averaging  the  the  same  equation 
of  motion  after  it  has  been  multiplied  by  the  appropriate  creation  or  destruction  operator 
and  then  scaling  the  resulting  equation  by  the  prefactor  in  Eqs.  (2.18-19) 


- 


<>iPh(r,t,X,+)  =tfpar-t(r)0ip/,(r ,  f,  A,  -f )  +  V'A*(r)0Op/l(rl  t) 


+  ”A'^A*(r)02p/»(r»M> +,*',-) 


(2.22) 


The  term  on  the  left  hand  side  of  the  equation  comes  from  pulling  the  time-dependent 
creation  operator  through  the  time  derivative.  The  zero-phonon-change  term  on  the  left 
hand  side  comes  the  term  in  the  last  sum  when  A  =  A';  for  this  term  the  two  explicit 
operators  (one  from  the  premultiplication  and  the  other  from  the  potential)  are  paired 
together  instead  of  with  two  unpaired  operators  in  the  amplitude-operator  and  the  ex¬ 
pectation  becomes  the  occupation  of  that  mode  times  the  zero-phonon-change  amplitude. 
The  correlated  contributions  are  again  negligble  because  of  unbalanced  factors  of  N~1. 


The  equations  of  motion  for  the  two-phonon-change  amplitudes  can  be  defined  in  a 
manner  similar  that  for  the  one-phonon-change  amplitudes,  leading  to  a  hierarchical  set 
of  equations  of  motion.  This  procedure  for  calculating  thermally  averaged  properties  is 
quite  general  for  any  system  in  which  a  few  degrees  of  freedom  are  coupled  to  an  infinite 
phonon  system.  The  approximations  we  have  made  up  to  this  point  have  been  made  just 
for  the  clarity  of  the  presentation.  For  instance,  the  inclusion  of  anharmonic  terms  in  the 
interaction  potential  just  complicates  the  equations  of  motion  for  the  thermally  averaged 
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amplitudes.  In  this  case  the  n-phonon-change  amplitudes  would  be  coupled  to  more  than 
just  the  (n  -  1)-  and  (n  +  l)-phonon-change  amplitudes.  Even  though  this  approach  is  an 
exact  solution  to  the  scattering  problem,  a  hierarchical  set  of  equations  of  motion  is  useful 
only  if  the  higher  order  phonon-change  terms  are  unimportant.  For  the  contributions  of 
the  n-phonon-change  terms  to  become  unimportant  in  the  limit  that  n  goes  to  infinity,  the 
coupling  between  the  phonons  and  the  scattering  particle  should  be  weak  as  is  the  case  for 
thermal  energy  helium  and  molecular  hydrogen  scattering. 


III.  ONE-PHONON  APPROXIMATION 

For  systems  in  which  the  inelastic  scattering  is  weak,  we  truncate  the  set  of  thermally 
averaged  amplitudes  to  include  only  zero-  and  one-phonon  changes.  The  equations  that 
describe  this  approximation  come  from  those  of  the  previous  section  with  the  higher  order 
phonon-change  amplitudes  set  to  zero.  The  equation  of  motion  for  the  zero-phonon-change 
amplitude  is  unchanged  because  the  interaction  potential  we  used  included  only  the  term 
linear  in  the  phonon  coordinates,  i.e. 


We  have  made  the  one-phonon-increase  and  one-phonon-decrease  amplitudes  look  similar 
by  defining  a  notation;  a  -  +  or  which  refer  to  a  net  phonon  increase  or  decrease, 
respectively,  after  the  interaction,  so  that  riy  +  =  n*  +  I,  =  n\  and  V^+  =  V*_  = 

In  this  notation  the  equations  of  motion  for  the  one-phonon-change  amplitudes  are  simple 
because  they  are  only  coupled  to  the  zero-phonon-changc  amplitude 

(‘^  ~aWA  “  HPart)  X'°)  =  lA*a(r)^0Pfc(r.O-  (3-2) 

In  the  approximation,  the  thermally  averaged  density  matrix  contains  contributions  from 
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only  the  zero-  and  one-phonon-change  amplitudes 

**(r,r\0  =  ^oph{r>t)rl}0ph{r' ,t)  +  ~  ^(r,*, \,o)il>lph{r',  t,  X,e)nXa.  (3.3) 

A  ,a 

If  the  density  matrix  is  evaluated  at  equal  postions  and  integrated  over  all  space  it  is 
constant  in  time 

j  d?r  P(r,r,t)  =  0.  (3.4) 

This  can  be  seen  by  taking  the  time  derivative  inside  the  spatial  integration  and  using 
the  equations  of  motion  to  replace  the  time  derivatives,  then  the  Hermitian  properties  of 
the  Hamiltonian  can  be  used  inside  the  integration  to  cancel  all  of  the  terms.  The  equal 
position  density  matrix  is  just  the  probability  density  for  the  scattering  particle  summed 
over  all  of  the  possible  states  of  the  lattice.  That  it  is  constant  in  time  means  that  the 
one-phonon-change  approximation  is  a  unitary  approximation.  The  same  arguments  apply 
to  any  finite  phonon-change  truncation  of  the  set  of  amplitudes  which  will  also  be  a  unitary 
approximation. 


IV.  STATIONARY  STATE  SCATTERING 

In  this  section  we  present  the  stationary  state  scattering  results  for  the  scattering 
probabilities,  discuss  the  problems  that  arise  in  stationary  state  scattering  theory  when 
the  zeroth  order  Hamiltonian  has  bound  states,  and  show  how  the  scattering  probabilities 
can  be  calculated  form  the  asymptotic  forms  of  the  scattering  state  wavefunctions.  Central 
to  the  derivation  of  these  results  (see  Appendix  B)  is  a  limiting  process  in  which  the  initial 
momemtum  uncertainty  of  the  wavepacket  goes  to  zero.  Taking  this  the  limit  allows  the 
scattering  to  be  calculated  in  terms  of  stationary  state  eigenfunctions  of  the  Hamiltonian. 

The  scattering  from  the  surface  can  be  described  in  terms  of  the  scattering  state 
solutions  of  a  Tippmann-Schwinger  equation  containing  the  static  surface  potential.  The 


state  xl  1  1S  a  scattering  state  with  outgoing  boundary  conditions  and  is  defined  using  a 
Lippmann-Schwinger  equation  with  the  advanced  free  particle  Green  function  Gq 

Xl!(r,k)  =  e'kr  +  j  dV  GS(r,rUk)V(rV_,(r>k)-  (4.1) 

Scattering  states  with  incoming  boundary  conditions,  x,  +  \  can  be  defined  in  the  same 
manner  using  the  retarded  Green  function.  The  bound  states  of  the  static  surface  potential 
are  defined  by 

(iW-£n,K)x(r,n,K)=0; 

j  d\  jX(r,n,K)|2  =  A.  (4.2) 

where  A  is  the  area  of  the  surface  and  n  indexes  the  bound  states  at  each  value  of  K.  A 
complete  set  of  states  consists  of  both  the  bound  states  and  either  the  incoming  or  the 
outgoing  scattering  states. 

The  propagation  of  a  particle  in  the  presence  of  the  static  surface  is  described  by  a 
retarded  static  surface  Green  function  Gj  defined  by 

E-Hpart  +  i*  Gr1(r,r',E)=6W(r~r'),  (4-3) 

together  with  the  boundary  conditions  that  in  the  limit  that  as  z  goes  to  infinity,  with  z' 
finite,  the  Green  function  behaves  like  an  outgoing  plane  wave.  We  use  this  Green  function 
to  write  scattering  states  in  the  presence  of  inelastic  scattering. 

The  scattering  states  in  the  presence  of  the  coupling  to  the  phonons  are  defined  by  a 
Lippmann-Schinger  equation  similar  to  those  with  no  coupling  by  using  an  outgoing  static 
surface  scattering  state  and  the  retarded  static  surface  Green  function.  The  Lippmann- 
Schwinger  equation  reduces  to  a  set  of  equations  coupling  the  different  multi-phonon  am¬ 
plitudes 


^OpUr’k)  =  Xl*'(r-k)  -  J  dV  G\(r,r',Ek)~  ^  ^ (r' ) ^ (r' ,  k.  A, cr),  (4.4) 
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(4.5) 


♦A14'1 

ph 


(r,  k,  A,  a)  =  -  J  d3r'  G\{r,r',Ek  -  aux)V*a( ,k). 


There  is  no  boundary  term  on  the  right  hand  side  of  this  last  equation  because  the  boundary 
conditions  are  chosen  so  that  the  incident  amplitude  in  the  one-phonon-change  states  is 
zero. 


A.  TRAPPING. 

The  bound  states  of  the  static  surface  potential  appear  as  poles  as  a  function  of  energy 
in  the  static  surface  Green  function.  Hence  the  one-phonon-change  amplitudes  diverge  as 
a  function  of  the  phonon  energy  as  the  final  state  energy  Ek  —  ou>\  becomes  equal  to  a 
bound  state  energy  with  parallel  wavevector  equal  to  K  —  crQ^.  These  divergences,  which 
arise  for  finite  temperature  in  any  finite  phonon-change  truncation  of  the  scattering  states 
and  for  zero  surface  temperature  in  any  approach,  are  due  to  trapping  on  the  surface.  A 
stationary  state  implies  a  constant  incident  flux  of  particles  which  leads  to  the  paradox 
that  the  number  of  trapped  particles  increases  continuously.  This  continual  build-up  of 
particles  means  that  there  is  no  time  independent  scattering  state  unless  there  is  absorption 
of  probabilty  from  the  frcm  the  bound  state  amplitudes;  this  absorption  (and  the  resolution 
of  this  pardox)  is  provided  by  the  imaginary  part  in  the  Green  function.  In  the  limit  that 
the  imaginary  part  of  the  Green  function  goes  to  zero  all  of  the  absorption,  which  is 
proportional  to  rj,  takes  place  from  the  bound  state  amplitudes  because  they  are  diverging 
proportional  to  tj~1.  The  limit  r?  — *  0  is  well  defined  and  produces  a  stable  limiting  value 
for  the  trapping  probability  (see  Appendix  B).  Even  though  these  scattering  states  involve 
a  loss  of  probability  due  to  trapping,  the  one-phonon  approximation  we  are  describing 
remains  unitary. 

The  stationary  state  scattering  wavefunctions  also  obey  a  coupled  set  of  Schrodinger 
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equations  with  the  infinitesimal  imaginary  part  included 

Hpart  -  Ey~  k)  =  ^  (4.6) 

L  A,cr 

Hpart  -  [Ev  -oux)  -  IT)  t^Jr.k ,A,o)  =  (4.7) 

Solving  these  equations  with  a  finite  value  of  rj  remedies  the  divergence  in  the  one-phonon- 
change  amplitudes  (see  Appendix  B).  In  practice  it  is  simple  to  find  a  finite  imaginary 
part  that  is  small  enough  that  that  the  solutions  of  these  equations  for  finite  rj  give  the 
same  results  as  those  in  the  limit  that  rj  goes  to  zero.  In  the  numerical  calculation  that 
we  describe  in  section  VI  we  choose  to  solve  the  Schrodinger  equations  instead  of  the 
Lippmann-Schwinger  equations  because  using  the  latter  would  require  either  storing  or 
recalculating  the  static  surface  Green  function,  G j. 

B.  SCATTERING  PROBABILITIES. 

To  describe  the  scattering  probabilities  in  the  presence  of  the  static  surface  potential, 
which  can  not  be  treated  as  a  perturbation,  we  have  to  calculate  the  matrix  elements  of 
the  wavepacket  with  the  outgoing  scattering  states  of  the  stat!  surface  potential.  This 
approach  is  quite  similar  to  that  of  the  distorted  wave  Born  approximation,  and  in  fact 
it  is  a  self-consistent  improvement  of  that  approximation.  In  particular,  the  4  robability 
to  scatter  into  a  state  with  a  wavevector  k j  is  given  by  the  square  of  the  matrix  element 
between  the  wavepacket  and  the  outgoing  scattering  state  solution  of  the  static  surface 

P( k/.O  =  ((*(0|x‘-,(r,k/))  (x'-V.k,)  *(0))^ 

=  \f  d\  X,',*(r,k/)^(r.01  (4  S) 

-f  ~N  YsnAj  dZr  X,'l*(r’k/)v A, o)' 
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In  Appendix  B  we  show  that  the  first  term  gives  the  elastic  scattering  and  the  remaining 
terms  give  the  inelastic  scattering  probabilities.  The  elastic  and  inelastic  scattering  in 
Eq.  (4.8)  can  be  written  in  terms  of  matrix  elements  of  the  the  interaction  potential 
with  respect  to  the  zero-phonon-change  scattering  state  amplitude  and  the  outgoing  static 
surface  scattering  states.  In  this  form  the  inelastic  amplitudes  look  like  a  typical  golden 
rule  rate  in  which  the  interaction  potential  causes  a  transition  from  the  initial  state  to 
the  final  state.  Note  that  this  approach  uses  the  incoming  zero-phonon-change  amplitude 
instead  of  the  incoming  static  surface  scattering  state. 


After  sufficiently  long  time  such  that  all  of  the  non-trapped  particles  have  left  the 
surface  the  probability  to  have  scattered  inelastically  into  a  final  state  with  wavevector  ky 
is  given  by  (see  the  discussion  following  (B6)) 


2 

(-1.9) 


n«,(k,)  /  i3’  rw/ifii'iCi'W 

zt  X,a 

(2tt)<5  (£k|  -  oux  -  £k/)  . 

This  result  is  identical  to  the  distorted  wave  Born  approximation  result  with  the  zero- 
phonon-change  amplitude  replacing  the  incident  static  surface  potential  scattering  state, 
i.e.  using  the  left  hand  side  of  Eq.  (4.4)  to  replace  the  first  term  on  the  right  hand  side 
of  Eq.  (4.4)  in  Eq.  (4.9).  This  expression  illustrates  why  this  approach  is  a  unitary 
extension  of  the  distorted  wave  Born  approximation  because  the  elastic  wavefunction  from 
which  the  inelastic  transitions  are  made  is  calculated  sclf-consistently  with  the  inelastic 
states  to  which  the  transitions  are  made.  Furthermore  we  see  why  the  distorted  wave 
Born  approximation  is  not  unitary  and  this  approximation  is;  in  the  distorted  wave  Born 
approximation  when  a  particle  scatters  inelastically  the  flux  stays  in  the  elastic  channel 
arid  continues  to  scatter  leading  to  an  overcounting  of  the  inelastic  scattering,  whereas  in 
this  approach  the  tlux  is  removed  from  the  elastic  channel  leading  to  a  flux  conserving 
approximation.  The  trapping  probabilities  are  the  probabilities  to  be  in  one  of  the  bound 


states  of  the  static  surface  potential  and  are  given  by  a  similar  expression  with  the  outgoing 
scattering  state  replaced  by  a  bound  state 


Pxnel{n,Kf)  /  d*r  Xl",*(r,n,K/)V^;(7(r)v»0p{l(r,kt) 


A, <7 


{2k)6  (£k<  -  awA  -  £niK/)  • 


(4.10) 


The  first  term  in  Eq.  (4.8)  gives  the  elastic  scattering  probability  which  can  be  shown 
to  be 


P«l(k/)  =  E  -  £k,K2*)^(2)(K,  -  K,  -  G) 


G 


k/.k.)  +  ^J  d3'  /iV  X|-|*(r,k/)E(r,r',£k,)^(r',k,-) 


(4.11) 


where  R( k/,kj  describes  the  amplitude  to  make  a  transition  from  the  incident  state  kt 
to  the  final  state  by  scattering  from  the  static  surface  potential.  It  is  proportional 
to  the  matrix  element  between  the  incoming  scattering  state  and  the  outgoing  scattering 
state.  The  sum  over  surface  reciprocal  lattice  vectors  G  makes  explicit  the  possibility  of 
elastic  diffraction  (or  inelastic  scattering  due  to  internal  excitations  of  the  molecule  as  we 
discuss  below)  in  the  elastic  scattering  probability.  The  second  term  in  the  absolute  value 
is  the  effect  of  the  inelastic  scattering  on  the  elastic  scattering;  it  is  written  in  terms  of  the 
self-energy  given  by 


E(r,r',  £k)  =  n\ov\Ar)Gri  (r. r',  Ek  -  <^A)TAcI(r').  (4.12) 

A, (7 

This  is  the  Born  approximation  to  the  self-energy  of  the  full  scattering  system  and  can 
be  understood  in  terms  of  its  constituents:  the  potential  causes  a  transition  to  a  one- 
phonon-change  state,  the  molecule  propagates  with  the  energy  Ek  -  owA  until  the  potential 
V^a  causes  a  transition  back  to  the  zero-phonon-change  state.  This  process  is  summed  over 
the  possible  phonons,  weighted  by  the  thermal  occupation  of  each  phonon  mode. 


The  inelastic  scattering  probabilities  can  also  be  directly  calculated  from  the  asymp¬ 
totic  form  of  the  one-phonon-change  scattering  states.  The  probabilty  to  scatter  into  a 
particular  final  state  is  given  by  the  flux  density  in  the  one-phonon-change  amplitude  far 
from  the  surface  divided  by  the  incident  flux  and  summed  over  all  phonon  modes  that 
satisfy  energy  and  momentum  conservation 

|2 


G  2t  X,<t 


4.13) 


—  (2* K(£k,  -  -  £k  )(2,r)  W(Kf-  -  aQA  -  G  -  Kf). 

The  trapping  probability  can  be  calculated  from  the  behavior  of  the  one-phonon-change 
amplitudes  in  the  limit  that  the  small  imaginary  part  in  the  Green  function  (or  Schrodinger 
equation)  goes  to  zero.  Calculating  the  scattering  state  solution  with  a  finite  imaginary 
part  leads  to  the  probability  of  absorption  due  to  trapping.  The  trapping  probability  is 
then  equal  to  the  probability  of  absorption  divided  by  the  incident  flux  summed  over  all 
the  phonon  modes  that  satisfy  parallel  momentum  conservation 


Ptnel{n,Kf)  =  lira  \  & r  |^(r>  K ,  A,  a) 

V  G  X  <y  J 


(4-14) 


( 2 7r ) <5 ( £7k ,  -  oux  -  £ln,K/)(2^)2<5(2)(Kl  -  aQA  -  G  -  K ,). 

Similarly  to  the  inelastic  scattering  probabilities,  the  elastic  scattering  probability  (includ¬ 
ing  diffractive  and  internal  transitions)  can  be  shown  to  be  equal  to  the  outgoing  flux 
divided  by  the  incident  flux 

=  ,'i.m  E  r  K;hr'k'>  -  rf 

G  2‘  (4.15) 

^(2 *)6(Ekt  -  Ekf)(2n)26^(K,  -G-Kf). 

These  last  three  expressions  (4.13-15)  are  used  to  calculate  the  scattering  probabilities 
once  we  have  solved  for  the  scattering  states  using  Eqs.  (4.6-7).  When  integrated  over  all 
momenta  and  summed  over  all  bound  states,  these  probabilities  sum  to  one  because  this 
is  a  unitary  approximation. 
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V.  APPROXIMATIONS  ON  THE  POTENTIAL 


Although  the  equations  presented  in  the  last  section  could  be  solved  for  an  arbitrary 
form  of  the  potential,  we  show  in  this  section,  starting  from  a  completetely  general  form  of 
the  interaction  potential,  what  assumptions  lead  to  a  Hamiltonian  for  which  the  scattering 
can  be  easily  calculated  and  what  form  the  equations  of  the  last  section  take  when  applied 
to  this  form  of  the  Hamiltonian.  In  particular  we  choose  a  Hamiltonian  for  a  flat  surface 
so  as  to  decouple  the  motion  parallel  to  the  surface  from  the  motion  perpendicular  to  it. 
This  decoupling  allows  a  much  simpler  numerical  calculation  of  the  wavefunctions. 


The  Hamiltonian  consists  of  the  kinetic  energy  of  the  scattering  particle,  the  lattice 
Hamitonian,  and  the  potential  that  couples  the  scattering  particle  and  the  lattice.  Later 
we  regroup  the  terms  in  the  pattern  of  Eq.  (2.3),  but  here  we  keep  all  of  the  coupling 
terms  together.  The  Hamiltonian  for  the  general  form  is 


H  =  Tpart  +  Hiat  +  V.  (5.1) 

The  kinetic  energy  of  the  particle,  Tpart ,  here  a  molecule,  consists  of  the  translational 
kinteic  energy  of  the  particle  plus  its  rotational  kinetic  energy.  We  neglect  vibrational  and 
electronic  degrees  of  the  molecule  because  the  excitation  energies  are  much  higher  than 


the  other  energies  in  the  scattering  problem.  Accordingly 

„  p2  L2 

Tpart  ~  2m  +  21  ’ 


(5.2) 


where  L  is  the  angular  momentum  operator  for  the  molecular  rotations  (a  three  dimensional 
vector  in  spite  of  its  being  upper  case)  and  /  is  the  moment  of  inertia.  The  orientation  of  the 
molecule  is  specified  by  the  angle  of  the  molecular  axis  with  respect  to  the  surface  normal 
0  and  the  angle  of  the  orientation  around  the  surface  normal  d>.  The  lattice  Hamiltonian 
is  the  same  as  it  was  in  the  previous  section  Eq.  (2.4) 

Hlat  =  '}2UJXaXa\- 
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The  displacements  of  the  lattice  atoms  can  be  expanded  in  terms  of  the  normal  modes  of 
the  lattice23 

u”  =  7^  E  "’«(*•  + «») '  =  -«»•  (s-i) 

Here  A  indexes  the  normal  modes,  is  the  wavevector  of  the  mode  parallel  to  the  surface, 
u>\  is  the  frequency  of  the  mode,  and  e(A,.zn)  is  the  polarization  vector,  which  depends 
on  the  distance  of  the  atom  from  the  surface  because  there  is  no  translational  symme¬ 
try  normal  to  the  surface.  The  polarization  vectors  obey  the  following  orthonormality 
condition 

i^e*(A,sn).e(V,zn)e',R"^-Q>^  =  «5A)A,.  (5.5) 

n 

Finally  a \  and  are  the  creation  and  destruction  operators  for  the  phonon  mode. 


A.  INTERACTION  POTENTIAL. 


In  general  the  interaction  between  the  scattering  particle  and  the  lattice  depends  on  the 
position  of  all  the  atoms  in  the  lattice  and  all  the  molecular  coordinates.  Since  the  mean 
square  displacements  of  the  lattice  atoms  are  small  compared  to  the  characteristic  lengths 
of  the  potential  a  Taylor  series  expansion  of  the  potential  in  terms  of  the  displacements  of 
the  lattice  atoms  should  converge  rapidly 


V(r,d,<t>,  {  Uj } )  =V0(r,  6, 4>)  + 


dV 


n,a 


5un,Q  l{u,}=0 
d2V 


Un,a 


+  2  hki  du*>°dumJ  l{uJ}=oUn,QUm^ 


(5.6) 


We  truncate  this  expansion  at  the  term  linear  in  the  phonon  coordinates,  consistent  with 
the  one-phonon-change  approximation  we  have  made  previously. 

Now  we  regroup  the  terms  in  the  Hamiltonian  so  that  all  the  terms  that  are  independent 
of  the  phonon  coordinates  are  grouped  together.24  In  addition  we  make  part  of  the  flat 
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surface  approximation  -  we  ignore  any  dependence  on  the  postion  of  the  particle  parallel 
to  the  surface  in  the  first  term  in  the  potential  - 


Hpart  ~  Tpart  -f-  Vq 


(5.7) 


The  rest  of  the  fiat  surface  approximation  is  made  in  Eq.  (5.10). 


The  linear  term  in  the  coupling  between  the  phonons  and  the  scattering  particle  can 
be  written  in  a  manner  that  makes  it  easy  to  see  all  of  the  approximations  that  we  make 

k*  =  £  £  73  £ <,(Q>+G)  evwg,  q A, *,  1, *)«.(*, + ->)■ 

G  a,z„  Viv  A  '  A 

(5.8) 

We  have  expanded  the  coordinates  of  the  lattice  atoms  in  terms  of  the  normal  modes  of 
the  lattice  Eq.  (5.4),  and  expressed  the  dependence  of  the  potential  on  the  position  of  the 
particle  in  the  plane  of  the  surface  by  Fourier  transforming  the  potential  with  respect  to 
the  sum  of  the  parallel  wavevector  of  the  phonon  mode  and  a  reciprocal  lattice  vector.  The 
dependence  of  the  potential  on  the  phonon  wavevector  and  the  reciprocal  lattice  vector  is 
given  by  an  integral  over  the  surface  Wigner-Sietz  cell  of  the  potential  summed  over  the 
lattice  sites  in  the  surface 


Vla(G,QA,z, *.,»,*)  =  ~  f  SR 

a  JSWSC  ^  dun,a 


{u,}=0 


(5.9) 


This  is  still  a  completely  general  form  for  the  linear  coupling  term  in  the  potential. 


The  approximation  that  we  make  for  the  coupling  potential  is  given  by  restricting  the 
form  that  Vla  takes  to 

vla(G,Qx,z,zn,e,4>)  =  -^g.o^.o^Qa)  (-—£—)  ■  (5-10) 

This  expression  includes  the  rest  of  the  flat  surface  approximation;  we  keep  only  coupling 
to  the  phonon  motion  normal  to  the  surface  6a<z,  do  not  allowing  umklapp  processes  (i.e. 
G  =  0),  and  do  not  allow  coupling  of  any  other  component  of  the  motion  with  rotation  of 


the  molecule  around  the  surface  normal.  For  simplicity  we  also  assume  that  the  scattering 
particle  only  couples  to  the  phonons  through  the  projection  on  the  top  surface  layer  6Zn  0. 
Finally  we  assume  that  the  dependence  of  the  potential  on  the  parallel  wavevector  of  the 
phonon  that  is  being  coupled  to  and  the  dependence  on  the  height  of  the  particle  above 
the  surface  and  its  orientation  are  separable.  The  dependence  on  the  parallel  wavevector 
is  assumed  to  be  Gaussian25 

tf(QA)  =  e~(Ql/2Qc).  (5.H) 

This  last  assumption  greatly  simplifies  numerical  calculations  because  the  only  dependence 
of  the  phonon  amplitude  on  the  wavevector  of  the  phonon  comes  in  through  the  net  energy 
transfer  to  the  motion  normal  to  the  surface  and  an  overall  scale  factor.  This  simplification 
is  discussed  further  below. 

This  separablity  assumption  can  also  be  thought  of  as  a  local  height  approximation, i.e. 
that  the  effect  of  the  phonons  is  to  locally  shift  the  origin  of  the  potential  without  changing 
its  shape  in  the  z-direction.  If  this  potential  is  then  expanded  in  terms  of  the  phonon 
coordinates  the  coefficient  of  the  linear  term  is  proportional  to  the  partial  derivative  of  the 
uncoupled  term  with  respect  to  the  height  above  the  surface.  Then  the  Fourier  transform 
of  H( Q)  gives  the  effect  of  a  lattice  vibration  at  one  surface  point  on  the  local  height  of 
the  surface  at  another  surface  point  as  a  function  of  the  distance  between  the  two  points. 
The  local  height  approximation  should  be  a  good  approximation  close  to  the  surface  where 
the  effective  height  is  dominated  by  the  closest  atoms;  further  from  the  surface  where  the 
potential  is  determined  by  more  and  more  atoms  the  approximation  should  break  done. 
Fortunately  the  inelastic  scattering  is  dominated  by  the  region  of  the  potential  close  to 
the  turning  point  so  that  the  breakdown  of  this  approximation  should  not  be  important. 
Within  this  approximation  the  neglect  of  the  layers  below  the  top  is  equivalent  to  assuming 
that  only  the  top  layer  affects  the  local  in  the  potential.  Again,  this  neglect  should  be  a 
good  approximation  close  to  the  surface. 
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The  interaction  potential  can  now  be  written  in  a  form  similar  to  Eq.  (2.6) 


(5.12) 


The  dependence  on  the  height  above  the  surface  and  on  the  orientation  with  respect  to 
the  surface  is  independent  of  the  phonon  to  which  the  scattering  particle  is  coupled.  The 
dependence  on  the  position  parallel  to  the  surface  comes  in  through  a  factor  that  conserves 
the  total  wavevector  of  the  system  parallel  to  the  surface.  The  factor  contains  all  the 
information  about  the  strength  of  the  coupling  to  each  phonon  mode 


^  =  V5a?^(a’0)*(<w- 

Below  we  see  that  M\  is  related  to  a  weighted  phonon  density  of  states. 


(5.13) 


B.  SCATTERING  STATES. 


Because  of  the  form  of  the  potential  we  have  assumed  the  parallel  wavevector  and  the 
azimuthal  quantum  number  are  conserved  and  hence  the  the  form  of  the  static  surface 
scattering  states  simplifies 


X(->(r,<U,k  =  C»K/ 


(5.14) 


The  remaining  part  of  the  wavefunction  is  independent  of  the  incident  parallel  wavevector. 
The  zero-phonon-change  amplitude,  on  the  other  hand,  still  depends  on  the  full  incident 
wavevector  because  the  coupling  to  the  phonons  depends  on  the  parallel  wavevector 


(5.15) 


The  scattering  probabilities  depend  on  the  matrix  element  of  the  derivative  of  the  static 
surface  potential  between  these  two  scattering  states. 
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The  details  of  the  phonons  come  into  the  scattering  probabilities  through  a  weighted 
projected  density  of  states 

C(Q,w)  =  j;  ^2riXa\Mx\2[2n)6{uj  -  awA)(27r)26(2)(Q  -  oQx).  (5.16) 


This  density  of  states  is  projected  onto  the  surface  layer  by  the  ea(A,0)  factor  in  the 
definition  of  M\,  and  weighted  by  the  thermal  occupation  of  the  modes,  by  the  amplitude 
of  the  mode,  and  by  the  the  (local-height)  phonon-cut-off  function  //( Q).  The  weighted 
phonon  density  of  states  indicates  how  likely  the  particle  is  to  interact  with  a  phonon  at  a 
particular  frequency  and  parallel  wavevector.  It  is  also  worth  a  reminder  at  this  point  that 
C(Q,u)  includes  both  phonon  creation  and  anihilation  events  through  a.  The  inelastic 
scattering  probability  is  given  by  the  product  of  three  factors:  the  inverse  of  the  incident 
velocity,  the  square  of  the  matrix  element  of  the  incoming  zero-phonon-change  state  and 
the  outgoing  static  surface  state  with  the  derivative  of  the  static  surface  potential,  [Iproof 
reader,  please  leave  commma,  thanks]  and  the  weighted  phonon  density  of  states 


Pined* fJf)  =jp-  (x1  ]{z,0,kzf,lf,m)\  dV°gjd-  ^l)ph(2^<k t.J*,"*)) 


(5.17) 


C(Kt-Kf,Et-Ef). 


This  equation  is  Eq.  (4.9)  rewritten  in  the  form  appropriate  for  the  flat  surface  approxima¬ 
tion.  One  nice  feature  of  this  approximation  is  that,  given  the  matrix  element  for  one  final 
z-component  of  the  wavevector,  the  scattering  probability  for  any  wavevector  with  that 
z-component  can  be  calculated  without  recalculating  the  scattering  state  wavefunctions. 
Since  these  wavefunctions  are  calculated  by  numerically  solving  a  Schrodinger  equation, 
the  savings  in  computer  time  can  be  considerable. 

To  rewrite  Eqs.  (4.6-7)  in  the  form  appropriate  for  the  flat  surface  approximation  it 
is  useful  to  define  some  energies  and  wavevectors.  The  incident  energy  due  to  both  the 
energy  in  center  of  mass  motion  normal  to  the  surface  and  the  rotational  energy  is  given 
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The  change  in  energy  in  motion  normal  to  the  surface  and  the  rotational  energy  due  to 
exciting  a  phonon  of  wavevector  Q  and  frequency  u  is  given  by 


A  Eg  =  Ezx  -  Ezf  =  u>  + 


(K,  -  Q)2  Kf 


(5.19) 


We  are  going  to  expand  the  motion  of  the  molecule  in  terms  of  spherical  harmonics  with 
respect  to  its  rotation  so  it  is  useful  define  the  kinetic  energy  the  incident  molecule  would 
have  in  each  rotational  state  far  from  the  surface 


kk  =  F  + *) 

2m  zt  21 


(5.20) 


It  is  also  useful  to  define  the  same  quantities  for  the  molecules  that  have  scattered  inelas- 


tically 


khi&Ez)  _  bzi 


^  =  SL-  A Ez. 
2m  2m 


(5.21) 


These  wavevectors  are  used  in  the  Schrodinger  equations  that  the  scattering  states  obey. 

We  can  take  advantage  of  the  fact  that  the  one-phonon-change  amplitudes  only  depend 
on  the  energy  in  the  motion  normal  to  the  surface  and  in  the  rotational  motion  up  to  a 
scale  factor  and  calculate  many  of  the  amplitudes  at  the  same  time  by  scaling  the  one- 
phonon-change  amplitudes  by  M\  and  defining  an  amplitude  that  only  depends  on  the 
change  in  that  energy 

Here  A  Eg  is  evaluated  at  u  =  cruj\  and  Q  =  If  Eq.  (4.13)  for  the  inelastic  scattering 

probability  is  written  using  the  wavefunction  defined  above,  there  is  a  factor  IA/>|2  associ¬ 
ated  with  the  sum  over  phonon  modes.  If  we  integrate  the  inelastic  scattering  probability 


over  possible  parallel  wavevectors  we  can  define  a  new  phonon  density  of  states  that  is  a 
measure  of  how  strongly  phonons  lead  to  a  particular  change  in  the  energy 
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(5.23) 


We  use  Eq.  (4.13)  rewritten  using  this  phonon  density  of  states  to  calculate  the  inelastic 
scattering  probabilities. 


C.  COUPLED-CHANNELS. 

For  molecular  hydrogen  scattering  the  potential  only  depends  weakly  on  the  orientation 
of  the  molecule  and  the  rotational  energy  splittings  are  comparable  to  the  scale  of  the  inci¬ 
dent  energies  so  that  it  is  useful  to  formulate  the  scattering  problem  in  a  coupled-channels 
approach  by  expanding  the  potential  in  Legendre  polynomials  and  the  wavefunctions  in 
spherical  harmonics.  These  expansions  can  be  truncated  after  a  few  terms.  The  expansion 
of  the  potential  and  its  derivative  are  given  by 

VoM)  =^Vj(r)fi(co5tf),  (5.24) 

l 

=  £  V/Pdcose).  (5.25) 

Since  the  wavefunctions  are  expanded  in  spherical  harmonics,  we  need  the  spherical  har¬ 

monic  matrix  elements  of  the  potential 

Vll'i*)  -  ]T  Vt..{z)  f  dn  y^flPrtcosOlYt'nV,*).  (5.26) 

l"  J 

To  simplify  the  wavefunctions  we  suppress  all  of  the  initial  conditions  in  v,  liting  down  the 
wavefunctions.  To  suppress  the  azimuthal  dependence  of  the  wavefunction  we  define  a 
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spherical  harmonic  with  the  azimuthal  dependence  removed 


(5.27) 


Using  this  function  we  can  write  the  expansions  of  the  wavefunctions  as 


(5-28) 

l 

1>'+p'h(z,0,AEz)  =  ±Ez)Ylm{6).  (5.29) 

l 

We  solve  for  these  wavefunctions  numerically  and  use  the  solutions  to  calculate  the  ther¬ 
mally  averaged  scattering  probabilities. 


The  wavefunctions  we  have  defined  above  obey  a  coupled  set  of  differential  equations 
that  are  the  formulations  of  Eqs.  (4.5-7)  in  the  approximations  that  we  have  made  for  the 
potential 
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The  first  term  on  the  left  hand  side  of  both  equations  contains  the  kinetic  energy  operator 
minus  the  kinetic  energ-  each  channel  of  the  molecule  far  from  the  surface,  the  second 
term  contains  the  static  s..r'  *ential  that  reflects  the  molecule  away  from  the  surface 

and  that  couples  different  rot..  channels  within  the  same  phonon-change  amplitude. 

The  right  hand  side  of  the  first  eq-  ition  is  the  coupling  of  the  zero-phonon-change  ampli¬ 
tude  to  all  of  the  one-phonon-change  amplitudes,  and  the  right  hand  side  of  the  second 
equation  is  the  coupling  of  each  one-phonon-change  amplitude  to  the  zero-phonon-change 
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amplitude.  These  terms  are  asymmetric  because  of  the  way  in  which  the  one-  phonon- 
change  amplitudes  where  scaled.  We  solve  these  equations  with  an  finite  imaginary  part 
that  is  small  enough  that  the  results  do  not  change  if  it  is  made  smaller,  typically  it  must 
be  much  smaller  than  all  of  the  widths  of  the  resonances  in  the  problem. 

By  discretizing  the  energy  mesh  for  the  one-phonon-change  amplitudes  these  equations 
could  be  solved  by  a  large  coupled  channels  calculation,  but  in  situations  in  which  the 
sum  of  the  distorted  wave  Born  approximation  results  are  not  much  greater  than  one  it 
can  be  faster  to  solve  these  equations  iteratively.  For  the  first  iteration  we  set  the  one- 
phonon-change  amplitudes  to  zero  and  solve  for  the  zero-phonon-  change  amplitude;  this 
gives  the  scattering  from  the  static  surface.  Then  we  use  this  result  for  the  zero-phonon- 
change  amplitude  to  calculate  the  one-phonon-change  amplitudes;  this  gives  the  distorted 
wave  Born  approximation  results  for  the  inelastic  scattering  probabilities.  Then  these 
amplitudes  are  used  to  calculate  the  zero-phonon-change  amplitude  and  the  calculation 
procedes  iteratively. 

D.  SCATTERING  PROBABILITIES. 

The  boundary  conditions  that  the  amplitudes  must  satisfy  are  that  the  amplitudes 
decay  to  zero  into  the  surface 

jlim)05i(2|I)1^(*,U£I)=O.  (5.32) 

The  zero-phonon-change  amplitude  far  from  the  surface  consists  of  a  unit  amplitude  in¬ 
coming  plane  wave  in  the  incident  channel  and  outgoing  plane  waves  in  all  the  rotational 
channels  that  have  positive  kinetic  energy  far  from  the  surface 


Those  rotational  channels  that  do  not  have  sufficient  energy  must  also  decay  away  from 
the  surface.  The  allowed  channels  for  the  one-phonon-  change  amplitudes  are  all  outgoing 
plane  waves 

limvi';(:,U£:)^i(A£z)e,i|l(A£l)2  ;  k2zl(AEz)  >  0.  (5.3-1) 

The  forbidden  channels  must  also  decay  to  zero  away  from  the  surface.  Solving  the  coupled 
Schrodinger  equations  Eqs.  (5.30-31)  for  these  boundary  conditions  give  the  scattering 
probabilities. 

The  elastic  and  rotationally  inelastic  scattering  probabilities  are  given  by  the  outgoing 
flux  in  each  channel  divided  by  the  incident  flux  (see  Eq.  (4.11)) 

Pl  =  YL\Rl\2-  (5.35) 

Kzt 

The  inelastic  scattering  probabilities  are  given  by  the  ratio  of  fluxes  times  the  phonon 
density  of  the  states  due  to  the  scaling  of  the  amplitudes  (see  Eq.  (4.9)  and  Eq.  (5.22)) 

Pl(AEz )  =  k^~^-C{/lEz)\Rl(AEz)\l .  (5.36) 

kzi 

The  trapping  probabilites  are  due  to  the  absorption  caused  by  the  infinitesimal  imaginary 
part  in  the  Schrodinger  equation,  and  are  given  by  the  probability  density  in  the  bound 
states  times  the  absorption  rate  (see  Eq.  (4.10)) 

1 0',;’ (*,/,£,,  -  En)\  V2  C(Ezt  -  En).  (5.37) 
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Since  this  is  a  unitary  approxmation  these  probabilities  all  sum  to  one 

+  =  >  <5-3S> 
l  l  J  {  1  * 

In  the  next  section  we  calculate  these  probabilities  for  the  case  of  HD  scattering  from 


copper. 
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These  probabilities  have  been  integrated  over  all  possible  parallel  momenta  that  lead 
to  the  same  change  in  the  energy  in  the  motion  of  the  molecule  normal  to  the  surface. 
To  recover  the  full  distribution  over  all  of  the  individual  final  states  it  is  just  necessary 
to  multiply  the  probability  that  has  been  calculated  for  the  change  in  normal  energy 
appropriate  to  that  final  state  by  the  integrand  of  Eq.  (5.23).  This  intregrand  is  the 
product  of  the  phonon  density  of  states  and  the  delta  function  that  determines  the  final 
normal  component  of  the  energy  from  the  properties  of  the  phonons. 

The  flat  surface  approximation  discussed  in  this  section  is  useful  because  it  greatly 
reduces  the  number  of  final  states  that  have  to  be  integrated  over.  This  reduction  arises 
from  the  decoupling  of  the  motion  normal  to  the  surface  from  that  in  the  plane  of  the 
surface.  It  is  a  useful  approximation  for  studying  rotationally  inelastic  scattering  from 
uncorrugated  surfaces  because  of  this  simplification.  It  is  not  a  useful  approximation  for 
quantitatively  calculating  scattering  probabilities  for  corrugated  surfaces,  but  by  calculat¬ 
ing  the  interaction  of  inelastic  scattering  with  rotational  transitions  it  should  be  possible 
to  qualitatively  discuss  scattering  from  a  corrugated  surface. 


VI.  RESULTS 


To  demonstrate  the  method  presented  in  this  paper  we  use  it  to  calculate  scattering 
probabilities  for  HD  scattering  from  copper  and  compare  the  results  with  the  results  of 
the  same  calculation  using  a  distorted  wave  Born  approximation.  The  main  purposes  of 
this  section  are  to  demonstrate  that  it  is  possible  to  carry  out  the  calculations  that  we 
have  outlined  in  this  paper  and  to  show  how  selective  adsorption  resonances  and  inelastic 
scattering  can  affect  each  other.  More  extensive  calculations  are  presented  in  a  subsequent 
paper  along  with  the  details  of  the  potential,  the  phonon  spectrum,  and  the  numerical 
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techniques  that  we  use. 


The  differences  between  the  distorted  wave  Born  aproximation  and  this  self-consistent 
one-phonon  approximation  are  most  obvious  at  a  selective  adsorption  resonance.  At  the 
energies  m.ar  th>  resonace  we  have  chosen  to  study,  the  HD  molecule  can  scatter  from  the 
uncorrugated  static  copper  surface  into  either  an  1=0  or  an  1  =  1  rotational  state.  Above 
the  resonance  but  near  it  most  molecules  leave  the  surface  rotationally  excited,  in  an  1  =  1 
state.  The  1=2  rotational  state  plays  a  significant  role  in  the  rigid  surface  scattering  when 
at  this  selective  adsorption  resonance  the  molecule  can  make  a  virtual  transition  into  the 
1=2  rotational  state  and  the  second  lowest  bound  state  of  the  potential.  When  the  molecule 
can  make  this  transition  it  tends  to  spend  a  long  time  near  the  surface  in  this  rotationally 
excited  state.  The  scattering  probabilities  are  greatly  affected  by  this  resonance  with  the 
elastic,  1=0,  scattering  probability  increasing  to  one  near  the  center  of  the  resonance.  These 
scattering  probabilities  can  be  seen  in  the  top  panel  of  Fig.l,  the  solid  curve  is  the  elastic, 
1=0,  scattering  probability  from  a  static  surface,  and  the  dashed  curve  in  the  rotationally 
inelastic,  1=1,  scattering. 

The  shape  of  the  elastic  scattering  probability  as  a  function  of  incident  energy  is 
characteristic  of  a  Fano  resonance.  This  is  not  surprising  because  selective  adsorption 
resonaces,  Fano  resonces,  and  Feshbach  resonances  all  result  from  the  coupling  of  a  bound 
state  into  a  continuum  of  states,  and  all  have  similar  lineshapes. 

The  elastic  and  rotationally  inelastic  probabilities  sum  to  one  for  all  incident  energies 
in  the  static  surface  calculation;  the  distorted  wave  Born  approximation  for  the  inelastic 
scattering  probabilities  does  not  alter  the  static  surface  scattering  probabilities.  For  the 
distorted  wave  Born  approximation  to  be  valid  the  inelastic  scattering  probability  should 
be  small  compared  to  one.  Far  from  the  resonances  the  inelastic  scattering  probability, 
seen  in  the  dotted  curve  in  the  top  pane!  of  Fig.l,  calculated  using  the  distorted  wave  Born 
approximation  is  about  as  large  as  it  can  be  for  this  approximation  to  remain  valid.  At 
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the  resonance  this  condition  is  clearly  not  satisfied  as  the  inelastic  scattering  probability 
excedes  unity.  The  total  probabilities  calculated  using  these  methods  range  from  about 
1.3  to  2.1.  The  width  of  the  Lorentzian  that  can  be  fit  to  the  enhancement  of  the  inelastic 
scattering  peak  is  the  same  as  the  width  of  the  Fano  resonace  that  can  be  fit  to  the  static 
surface  scattering  probabilities. 

The  bottom  panel  of  Fig.  1  for  comparison  shows  the  results  of  the  same  calculation 
done  using  the  method  developed  in  this  paper.  Since  this  approximation  is  unitary  the 
scattering  probability  is  one  for  all  incident  energies;  numerically  unitarity  is  satisfied  to 
the  same  accuracy  as  it  is  for  the  static  surface  calculation  of  the  elastic  and  rotation- 
ally  inelastic  scattering  probabilities.  In  this  approximation  the  elastic  and  rotationally 
inelastic  scattering  probabilities  are  affected  by  the  inelastic  scattering  probabilities.  The 
peak  in  the  elastic  scattering  probability  is  strongly  reduced,  much  more  strongly  than 
is  the  dip  in  the  rotationally  inelastic  scattering  probability.  This  suggests  that  dips  due 
to  selective  adsorption  resonances  should  be  easier  to  observe  experimentally  than  peaks. 
The  width  of  the  enhancement  peak  in  the  inelastic  scattering  has  increased  due  to  the 
interaction  between  the  selective  adsorption  resonace  and  the  inelastic  scattering;  the  in¬ 
crease  is  roughly  a  factor  of  two.  This  increase  in  the  resonance  width  indicates  that  the 
inelastic  lifetime  of  the  resonace  is  comparable  to  the  inherent  rotational  width. 

VII.  SUMMARY 

The  central  features  discussed  in  this  paper  are:  (1)  when  the  scattering  is  weak, 
we  justify  expanding  the  scattering  wavefuctions  in  terms  of  n-phonon-change  operator- 
amplitudes,  (2)  these  n-phonon-change  operators  amplitudes  are  thermally  averaged  to 
describe  the  averaged  scattering  ior.  terms  of  thermally  averaged  amplitudes,  (3)  the  srat- 
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tering  probabilities  are  calculated  using  stationary  scattering  states  when  the  interaction 
time  is  long  compared  to  the  characteristic  times  in  the  substrate,  (4)  the  calculation  can 
be  simplified  for  some  specific  assumptions  about  the  interaction  potential,  (5)  and  the 
inelastic  scattering  and  selective  adsorption  resonances  affect  each  other.  The  first  three 
points  reflect  what  we  feel  is  the  correct  way  to  treat  thermal  energy  helium  and  molecular 
hydrogen  scattering  from  surfaces  on  which  they  physisorb.  The  fourth  point  details  the 
simplication  of  the  scattering  calculation  for  a  specific  system  in  which  the  corrugation  of 
the  surface  does  not  play  an  important  role  and  the  fifth  point  is  a  result  of  applying  these 
approximations  to  molecular  hydrogen  scattering. 

1.  When  the  inelastic  scattering  is  sufficiently  weak  that  in  a  typical  scattering  event 
only  a  few  phonons  are  created  or  destroyed,  the  scattering  process  can  be  reasonably 
be  described  in  terms  of  how  the  occupation  of  the  lattice  has  changed.  Helium  and 
molecular  hydrogen  scattering  at  thermal  energies  from  surfaces  on  which  they  physisorb 
satisfy  the  weak  inelastic-scattering  criterion  due  to  their  low  mass  and  their  weak  in¬ 
teraction  with  the  surface.  Since  the  elastic  scattering  is  a  quantum  mechanical  process 
and  is  both  observable  and  distinguishable  from  the  inelastic  scattering,  it  is  necessary  to 
treat  the  elastic  scattering  quantum  mechanically  and  seperately  from  the  inelastic  scat¬ 
tering.  It  is  also  important  to  treat  the  inelastic  scattering  in  terms  of  the  changes  in  the 
phonon  modes  because  the  inelastic  scattering  probabilities  reflect  the  discrete  nature  of 
the  phonon  excitations  of  each  mode.  Doing  the  calculation  self-consistently  allows  the 
range  of  validity  of  the  calculation  to  be  extended  to  resonant  elastic  scattering  situations 
in  which  a  perturbative  approach  would  break  down. 

2.  In  calculating  scattering  probabilities  it  is  useful  to  be  able  to  thermally  average 
the  results  without  having  to  do  an  explicit  ensemble  average  of  calculated  scattering 
probabilities.  In  this  paper  we  have  presented  a  method  of  doing  the  thermal  averaging 
by  calculating  the  scattering  in  terms  of  thermally  averaged  n-phonon-change  amplitudes. 
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This  part  of  the  calculation  is  exact  and  is  useful  in  any  situation  in  which  the  expansion 
in  terms  of  phonon  exchanges,  discussed  in  the  previous  section,  is  applicable. 

3.  The  scattering  probabilities  for  these  systems  can  be  calculated  using  stationary 
state  scattering  wavefunctions  because  the  the  resolution  with  which  the  energy  of  the 
incident  particle  is  defined  is  much  smaller  than  typical  phonon  frequencies  in  the  lat¬ 
tice.  Stationary  state  calculations  are  useful  both  for  their  simplicity  and  because  they 
emphasize  the  quantum  mechanical  nature  of  the  scattering  process  that  is  of  interest  in 
the  scattering  experiments.  Since  scattering  experiments  are  performed  for  a  better  un¬ 
derstanding  of  the  gas-surface  interaction  potential  and  quantum  mechanical  scattering, 
due  its  discrete  nature,  is  usually  more  sensitive  to  the  details  of  the  potential,  stationary 
state  calculations  are  better  suited  to  discriminate  between  possible  potentials. 

4.  The  flat  surface  approximation  is  used  because  it  speeds  up  the  calculation  while  still 
including  the  some  of  the  important  aspects  of  the  potential  and  allowing  the  possibility  of 
studying  resonance  phenomena.  It  is  an  approximation  in  which  we  neglect  the  corrugation 
of  the  surface  while  still  treating  the  motion  of  the  particle  parallel  to  the  surface.  The 
simplification  comes  from  the  seperability  of  the  motion  parallel  and  perpendicular  to  the 
surface.  This  approximation  also  allows  separate  calculation  of  the  importance  of  the 
details  of  the  potential  and  of  the  importance  of  the  details  of  the  phonon  spectrum  on 
the  scattering  probabilities. 

5.  Our  preliminary  calculations  show  that  selective  adsorption  resonances  increase  the 
inelastic  scattering  at  resonance  condtions  and  that  inelastic  scattering  broadens  selec¬ 
tive  adsorption  resonces.  This  last  result  will  not  emerge  from  a  low  order  perturbative 
approach.  Future  work  will  be  directed  toward  studying  in  more  detail  how  inelastic  scat¬ 
tering  and  selective  adsorption  esonanccs  affect  each  other  as  well  as  how  finite  substrate 
temperature  affects  the  scattering  process. 
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APPENDIX  A:  THERMAL  AVERAGING  (DUE  TO  MDS) 

In  this  appendix  we  discuss  the  thermal  averaging  procedure  that  we  use  in  section 
II  of  this  paper.26  We  show  that  in  the  limit  that  the  number  of  modes  in  the  substrate 
goes  to  infinity  the  contribution  of  correlated  terms  to  the  density  matrix  goes  to  zero. 
When  the  correlated  contributions  are  unimportant  the  thermal  averages  of  products  of 
operators  can  be  written  as  products  of  thermal  averages  of  the  operators.  In  particular  we 
are  interested  in  the  thermal  averages  of  the  amplitude-operators  that  describe  the  time 
evolution  of  the  amplitude  of  the  scattering  particle;  these  are  the  thermally  averaged 
n-phonon-change  amplitudes  we  use  to  calculate  the  scattering  probabilities. 

The  correlated  contributions  to  the  thermal  expectation  value  are  unimportant  because 
they  are  inversely  proportional  to  the  size  of  the  lattice.  Correlations  give  rise  to  this  factor 
of  N~l  by  removing  a  sum  over  a  number  of  terms  proportional  to  the  size  of  the  system. 
The  form  of  the  interaction  potential,  Eq.  (2.12),  leads  to  both  a  sum  over  the  modes  of 
the  lattice  and  a  factor  of  N~ */2  associated  with  each  phonon  creation  and  anihilation 
operator  when  they  occur  in  the  density  matrix  Eq.(2.15). 

Supressing  the  spatially  dependent  functions  associated  with  each  phonon  operator 


leads  to  the  following  form  for  all  of  the  terms  in  the  left  hand  side  of  Eq.  (2.15) 


with  arbitray  numbers  of  creation  and  destruction  operators.  The  only  terms  in  this 
expression  that  are  non-zero  are  the  terms  which  have  no  unpaired  phonon  operators. 
When  all  the  operators  are  paired  there  are  the  same  number  of  sums  over  the  phonon 
modes,  as  there  are  factors  of  TV-1,  each  sum  having  37V  terms.  The  terms  that  contribute 
can  be  in  the  form 
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Here  all  of  the  phonon  operators  are  paired  in  the  same  way  that  they  are  in  the  right  hand 
side  of  Eq.  (2.15),  which  has  been  written  in  terms  of  operators  and  amplitude  operators 
that  have  the  pairing  already  explicit.  The  uncorrelated  terms  are  those  in  which  we 
replace  the  average  of  the  product  of  paired  operators  by  the  product  of  the  averages  of 
the  paired  operators 
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In  this  expression  all  factors  of  TV-1  are  balanced  by  sums  over  the  lattice  modes.  On  the 
other  hand,  correlated  terms  have  at  least  two  of  the  paired  operators  paired  with  each 
other  and  are  of  the  form 
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Since  each  of  the  correlated  terms  have  at  least  one  factor  of  TV-1  associated  with  it  that 
is  not  balanced  by  a  sum  over  the  phonon  modes  these  terms  do  not  contribute  in  the  limit 
that  the  number  of  lattice  atoms  goes  to  infinity.  The  vanishing  of  the  correlated  terms 
means  that  we  can  write  Eq.  (2.15)  as  Eq.  (2.16),  and  then  define  the  n-phonon-change 
amplitudes  that  we  use  to  calculate  the  thermally  averaged  scattering. 


As  an  example  of  how  terms  from  the  right  hand  side  of  Eq.  (2.15)  contribute  to  the 
right  hand  side  of  Eq.(2.l6)  let  us  consider  one  of  the  terms  from  the  square  of  the  zero- 


phonon-change  amplitude-operator  (see  Eq.  (2.14)  for  the  definition  of  these  operators) 
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First  consider  the  uncorrelated  contributions,  those  in  which  Aj  /  A2,  rewriting  the  expo¬ 
nential  of  a  sum  as  a  product  of  exponentials 

=  jf  2Z  tv  r!  !  n-e-'^ (A6) 

Ai  A2^Aj  ^in>}  Ut 

Now  write  the  sum  over  all  sets  of  occupation  numbers  of  the  products  over  all  modes  as 
the  product  over  all  modes  of  the  sums  of  the  occupations  for  each  mode 


rn=nL 

{n,}  1  t  n,  =0 


(417) 


All  the  factors  in  the  numerator  other  than  these  containing  Aj  and  A2  cancel  similar 
factors  in  the  partition  function  leaving 

x  x  V'°°  ^  ne~PUJxin  V°°  ^  ne~^x2n 
r.Aj *  ,A2  E,n  = 0ne  1  L^n- 0ne 
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These  two  ratios  just  give  the  thermal  expectation  of  the  occupation  of  each  of  the  modes 
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These  and  all  the  other  uncorrelated  terms  terms  are  the  contributions  to  the  right  hand 
side  of  Eq.  (2.16).  Now  consider  the  correlated  terms;  the  calculation  for  these  terms 
remains  the  same  until 
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Although  this  result  is  of  a  different  form  than  the  terms  in  Eq  (A8),  there  is  an  additional 
N~l  factor  associated  with  losing  a  sum  over  the  phonon  modes  so  that  the  difference  is 
unimportant.  This  factor  causes  the  contributions  of  the  correlated  terms  to  go  to  zero  in 
the  limit  that  the  number  of  lattice  atoms  goes  to  infinity.  The  restricted  sum  in  Eq.  (A9) 
can  be  made  unrestricted  in  the  same  limit  because  the  error  is  also  proportional  to  TV-1. 

APPENDIX  B:  STATIONARY  STATE  SCATTERING 

In  this  appendix  we  outline  the  derivation  of  the  results  in  section  IV  of  this  paper. 
First  we  show  how  to  calculate  the  scattering  probabilities  in  terms  of  matrix  elements 
of  the  zero-phonon-change  scattering  state  by  taking  the  limit  that  the  initial  momentum 
uncertainty  of  the  wavepacket  goes  to  zero.  Then  we  show  that  the  same  scattering  prob¬ 
abilities  can  be  calculated  from  the  asymptotic  form  of  the  zero-  and  one-phonon-change 
scattering  states. 

The  derivation  of  the  results  in  this  appendix  requires  four  limits  to  be  taken  in  this 
order:  the  imaginary  part  rj  that  determines  the  boundary  conditions  goes  to  zero,  the 
wavepacket  start  far  enough  from  the  surface  that  it  is  not  interacting  with  it,  the  initial 
momentum  uncertainty  of  the  wavepacket  goes  to  zero,  and  the  final  time  minus  the  initial 
time  (t f  —  <t)  goes  to  infinity.  The  first  two  and  the  last  limits  are  properties  of  description 
of  the  calculation;  the  imaginary  part  has  to  be  small  enough  that  the  normalization  of 
the  incident  plane  is  not  affected  by  it,  we  choose  to  have  the  particle  start  in  a  state  that 
is  independent  of  the  surface,  and  want  to  measure  the  state  of  the  particle  after  it  has 
ceased  interacting  with  the  surface. 

The  third  limit  is  a  property  of  the  scattering  system  and  may  not  always  be  valid.  To 
take  this  limit  (the  initial  momentum  uncertainty  going  to  zero)  the  interaction  time  of 


the  scattering  process  should  be  longer  than  the  characteristic  time  scales  of  the  phonons. 
The  interaction  time  is  set  by  the  time  difference  from  the  time  when  the  wavepacket  first 
interacts  with  the  phonons  to  the  time  when  all  of  the  incident  wavepacket  has  entered 
the  interaction  region.  A  lower  bound  on  the  interaction  time,  set  by  the  energy-time 
uncertainty  principle,  is  Planck’s  constant  divided  by  the  energy  resolution  of  the  incident 
wavepacket.  A  lower  bound  on  this  lower  bound  is  set  by  the  energy  resolution  of  the 
scattering  experiment.  Since  the  times  scales  of  the  phonons  are  their  oscillation  periods 
and  the  phonon  frequencies  go  continuosly  to  zero,  there  are  always  phonons  in  the  surface 
with  time  scales  longer  than  the  interaction  time.  For  this  reason  stationary  state  scattering 
is  not  be  able  to  describe  all  of  the  scattering  processes,  but  if  the  interaction  time  is 
longer  than  most  relevant  time  scales  and  if  the  low  frequency  phonons  do  not  dominate 
the  scattering  process  it  should  be  able  to  describe  the  scattering  within  the  resolution  of 
the  experiment.  In  particular,  if  low  energy  modes  are  not  important  in  the  stationary 
state  scattering  calculation  they  are  not  be  important  in  a  wavepacket  scattering  situation 
because  a  wavepacket  do  not  strongly  excite  modes  for  which  the  oscillation  period  is 
longer  than  the  interaction  time. 

To  derive  the  martix  elements  for  the  scattering  probabilities  we  expand  the  initial 
wavepacket  in  terms  of  the  scattering  state  solutions  in  the  one-phonon-change  approx¬ 
imation.  Then  we  use  the  spectral  representation  of  the  Green  function  to  convert  the 
time  dependent  exponentials  into  energy-conserving  delta  functions  in  the  limit  that  the 
time  difference  goes  to  infinity.  Finally  the  amplitude  factors  that  describe  the  initial 
wavepacket  become  a  delta  function. 

The  initial  wavepacket  can  be  expanded  in  incoming  plane  waves 

rf(r,<)  =  /  (0>  (Bl) 

where  the  time  dependence  gives  the  free  motion  of  the  particle.  The  expansion  coefficient, 
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a(k),  is  sharply  peaked  around  the  incident  wavevector  kt.  The  boundary  conditions  that 
the  Lippinann-Schwinger  equations  for  the  zero-  and  one-phonon-change  amplitudes  obey 
allow  us  to  expand  the  initial  wavepacket  in  terms  of  these  scattering  states  with  the  same 
coefficients  as  the  expansion  in  plane  waves 


=  J  a(k)V^(r,k)e  t£k(£  £'>,  (B2) 

0lPh(r,<,A,CT)  =  J  a(k)0(1pji(r,k,A,a)e-l£:^£_£>).  (£3) 

This  time  dependence  gives  the  full  time  dependence  of  the  scattering  particle  in  the 
one-phonon-change  approximation.  Inserting  the  equation  for  the  one-phonon-change  am¬ 
plitude  into  Eq.  (4.7),  gives  the  expansion  of  the  state  of  the  scattering  particle  at  any 
time  in  terms  of  outgoing  scattering  states  of  the  static  surface.  To  write  the  scattering 
probability  in  the  form  of  Eq.  (4.9)  the  one-phonon-change  amplitudes  are  written  in 
terms  of  the  zero-phonon-amplitude  using  the  Lippmann-Schwinger  Eq.  (4.5),  and  the 
static  surface  Green  function  is  written  in  terms  of  its  spectral  representation 

rr(rr,  n  f  d^k  x<->((i,k))x(-'(r',k )*  f  Sk  ^x(r,n,K)x(r',K ,n)* 

l[  ’  ’ ]~1  (203  E-Ek  +  it,  +J  (2tt)2  L*  E  -  EnK  +  irj  ‘ 

(BA) 

The  sum  over  n  in  the  second  sum  in  the  spectral  representation  is  restricted  to  run 
only  over  the  bound  states  for  a  particular  value  of  K;  since  all  of  the  Enj&  must  be 
less  than  zero,  some  values  of  K  have  no  bound  states  for  sufficiently  large  energy  in 
the  motion  parallel  to  the  surface.  The  normalization  of  the  scattering  states  allows  one 
three-dimensional  spatial  and  one  three  dimensional  k-space  integral  to  be  done  trivially, 
leaving 
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This  expression  is  valid  for  any  time,  in  particular  the  limit  that  tj  —  f,  — ►  oo  in  which  we 
calculate  the  scattering  probabilities. 


Expanding  the  square  of  the  absolute  value  as  a  product  of  complex  conjugates  gives 


e~i(Ek~i7iJx~EVf+trj)(tf-tl 

Ek  T^\  -  Ek,  +  IV 

el(Ev 

Ek'  -  ~  Ek,  ~  "1 


I  d3r  x{ 

I  d3r  x(''*(r,k/)V"A<J(r)*^pji(r,k/) 


Each  of  the  factors  with  time  exponentials  become  energy  delta  functions  in  the  limit  that 
the  time  difference  goes  to  infinity  and  the  imaginary  part  goes  to  zero,  as  can  be  seen  by 
converting  each  expression  into  a  time  integral.  Because  the  interaction  potential  conserves 
parallel  momentum,  we  can  factor  from  each  of  the  matrix  elements  a  parallel  momentum 
conserving  delta-function 

j(2n)26^(Kf  +oQx  +  G-K),  (B7) 

which  gives  unity  when  integrated  over.  The  expansion  coefficients  have  chosen  so  that  the 
overlap  integral  over  k  is  zero  if  one  of  the  arguments  is  displaced  by  a  surface  reciprocal 
lattice  vector 

f  d3k 

J  a(k)a(k  +  G)  =  0,  (58) 

and  their  modulus  squared  becomes  a  delta  function  as  the  initial  momentum  uncertainty 
goes  to  zero.  Integrating  over  the  delta  functions  gives  the  inelastic  scattering  probability 
result,  Eq.  (4.8).  Replacing  the  outgoing  static  surface  scattering  state  by  a  bound  state 
and  performing  the  same  manipulations  gives  Eq.  (4.9)  for  the  trapping  probability. 
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The  derivation  of  the  elastic  scattering  probability  proceeds  in  the  same  manner,  com¬ 
plicated  by  the  presence  of  outgoing  states  scattered  from  the  static  surface.  The  elastic 
scattering  from  the  static  surface  is  described  by  a  reflection  coefficient,  R( k^,k,)  that 
includes  the  possibility  of  diffraction 

/  i1rX'-|,(r,l/)*l'l(r,k,l  =  ^fl(k/lk,}(2»)3«(2)(K,-K/-G)^i(£,-E/).  (BO) 
J  G 

Inserting  Eq.  (4.5)  into  Eq.  (4.4)  gives  an  equation  for  the  zero-phonon-change  scattering 
state  in  terms  of  itself  and  the  self-energy,  Eq.  (4.11).  The  elastic  scattreing  probability 
is  given  by  inserting  the  resulting  expression  into  the  first  term  on  the  right  hand  side  of 
Eq.  (4.7).  After  factoring  out  delta  functions  from  all  of  the  terms  using  the  same  tricks 
that  we  used  for  the  inelastic  scattering  probabilities  for  the  terms  in  which  the  self-energy 
appears,  the  derivation  of  Eq.  (4.10)  procedes  in  much  the  same  manner  as  the  derivation 
of  Eq.  (4.8).  The  second  term,  containing  the  self-energy,  in  Eq.  (4.10)  cancels  the  part 
of  the  elastic  scattering  intensity  that  has  been  lost  to  inelastic  scattering.  The  energy 
conserving  delta  function  restricts  possible  outgoing  states  to  those  that  have  the  same 
energy  as  the  incident  particle,  i.e.  elastic  scattering. 

Rotational  degrees  of  freedom  only  complicates  these  derivations  by  adding  a  rotational 
state  subscript  to  almost  all  of  the  factors.  The  rotational  transitions  can  be  treated  in 
the  same  way  that  diffraction  was  treated  in  the  elastic  scattering  probability. 

The  rest  of  this  appendix  is  concerned  with  the  calculation  of  these  scattering  prob¬ 
abilities  from  the  asymptotic  forms  of  the  scattering  states.  For  the  inelastic  scattering 
probabilities  we  need  the  asymptotic  form  of  the  static  surface  Green  function  in  the  limit 
that  the  z-component  of  the  first  argument  goes  to  infinity 
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where  the  K/  integration  is  restricted  to  those  values  for  which  k'z 2  >  0.  Using  this  form 
and  Eq.  (4.5),  the  asymptotic  form  of  the  wavefunction  is 


(  +  )  i.  \  R„t 


f  d\'  x'-'-(r',k')Klo(rr^;k(r',k) 


K'  =  K  -aQA, 


k!z 2  =  2m(£],  —  oujx)  ~  Aw2 


{B 11) 


The  flux  density  in  this  state  for  large  distances  (squaring  it  and  multiplying  by  the  outgo¬ 
ing  velocity)  divided  by  the  incident  velocity  gives  the  contribution  from  each  mode  to  the 
inelastic  scattering  probability.  Summing  this  result  over  all  of  the  phonon  modes  weighted 
by  the  thermal  occupation  of  each  mode  and  multiplying  by  energy  and  parallel  momen¬ 
tum  conserving  delta  functions  gives  the  result  for  the  inelastic  scattering  distribution  Eq. 
(4.12).  A  similar  treatment  of  the  elastic  scattering  probability  yields  Eq.  (4.10).  These 
results  show  that  the  scattering  probabilities  can  be  given  by  the  flux  in  each  outgoing 
channel  of  the  zero-  and  one-phonon-change  amplitudes  divided  by  the  incident  flux. 


The  expression  for  the  trapping  probability,  Eq.  (4.9)  can  be  calculated  from  the 
behavior  of  the  static  surface  Green  function  in  the  limit  that  the  infinitestimal  imaginary 
part  goes  to  zero 

„r,  ,  ™  r  d2K'  x(r,n,K')x(r',n,K')*  /n,^ 

r;— ‘0,E— »£„,k  1  J  (2zr)2  E  -  £„,K'  +  "1 

Using  this  expression  the  one-phonon-change  amplitude  at  the  bound  state  energy  becomes 

1  f  d?  K' 

Jim  ^((r.k.k.a))  =  - J  ^  x(r.n.K') 

[  J  <iVx(r',n,Kr‘'A„(rr<;,i(r',k)  {Bn) 

K'  =K  -  oQx, 
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Since  the  bound  state  wavefunctions  are  normalized  to  the  area  of  the  surface,  Eq.  (4.2), 
squaring  the  amplitude,  integrating  the  result  over  all  spare,  dividing  by  the  area  of  the 
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surface,  and  multiplying  by  the  imaginary  part  squared  gives  the  contribution  to  the 
trapping  for  each  phonon  mode.  Summing  over  all  the  phonon  modes  weighted  by  the 
thermal  occupation  and  multiplting  by  energy  and  parallel  momentum  conserving  delta 
functions  gives  the  result  for  the  trapping  probability  Eq.  (4.13).  This  result  can  also  be 
derived  from  the  absorption  of  probability  from  the  system  due  to  the  small  imaginary  part 
and  trapping.  In  this  case  we  have  to  evauate  Eq.  (B13)  for  energies  near  the  bound  state, 
square  the  amplitudes,  multiply  by  the  absorption  rate,  rj,  and  integrate  over  both  all 
space  and  the  energies  near  the  bound  states.  For  small  enough  rj  all  the  matrix  elements 
are  constant  so  the  energy  integral  is  just  the  integral  of  a  Lorentzian,  which  gives  the 
additional  factor  of  t]  that  appears  in  Eq.  (B13). 

The  only  difficulty  that  arises  in  using  these  results  for  the  flat  surface  approximation 
is  that  the  the  spectral  representation  of  the  Green  function  no  longer  has  the  same  form 
as  it  does  in  Eq.  (B4)  because  the  motions  parallel  and  perpendicular  to  the  surface 
decouple.  In  particular  in  the  bound  state  terms  there  are  the  same  number  of  bound 
states  for  every  parallel  momentum,  some  of  which  have  a  positive  total  energy.  In  a  real 
system  the  slightest  coupling  of  the  two  motions  would  couple  these  positive  energy  states 
to  the  continuum  and  they  are  no  longer  bound  states.  This  defect  of  the  flat  surface 
approximation  is  not  important  for  one-phonon-change  scattering  at  normal  incidence  but 
would  complicate  the  interpretation  of  the  results  for  off-normal  scattering  conditions  when 
the  final  parallel  momentum  can  be  high. 
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FIGURE  CAPTION 


Fig.  1.  A  comparison  of  two  methods  for  calculating  scattering  probabilities.  The 
system  is  HD  scattering  at  normal  incidence  from  a  zero  temperature  copper  surface.  The 
top  panel  shows  as  a  function  of  incident  energy  a  static  surface  calculation  of  the  elastic 
(solid  curve,  /  =  0),  and  rotationally  inelastic  (dashed  curve,  /  =  1)  and  a  distorted  wave 
Born  approximation  calculation  of  the  total  phonon  inelastic  scattering  probability  (dotted 
curve).  The  bottom  panel  shows  the  same  probabilities  calculated  using  the  self-consistent 
one-phonon  approximation  outlined  in  this  paper.  The  structure  in  these  curves  is  due 
to  a  selective  adsorption  resonance  in  which  the  HD  molecule  makes  a  virtual  rotational 
transition  to  an  /  =  2  state  in  the  second  bound  state  of  the  static  surface  potential.  The 
three  probabilities  in  the  top  panel  sum  to  a  total  probability  greater  than  one,  especially 
at  the  resonance,  while  the  three  probabilities  in  the  bottom  panel  sum  to  one  at  all  values 
of  the  incident  energy. 


